MEPP2 Project
Interpolation.hpp
Go to the documentation of this file.
1 // Copyright (c) 2012-2019 University of Lyon and CNRS (France).
2 // All rights reserved.
3 //
4 // This file is part of MEPP2; you can redistribute it and/or modify
5 // it under the terms of the GNU Lesser General Public License as
6 // published by the Free Software Foundation; either version 3 of
7 // the License, or (at your option) any later version.
8 //
9 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
10 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
11 #pragma once
12 
13 #include <vector>
14 
15 namespace FEVV {
16 namespace Operators {
17 namespace Geometry {
18 
34 template< typename GeometryTraits >
35 static inline bool
36 is_in_tetrahedron(const typename GeometryTraits::Point &a,
37  const typename GeometryTraits::Point &b,
38  const typename GeometryTraits::Point &c,
39  const typename GeometryTraits::Point &d,
40  const typename GeometryTraits::Point &Pos,
41  const GeometryTraits &gt)
42 {
43  typedef typename GeometryTraits::Vector Vector;
44  typedef typename GeometryTraits::Scalar Scalar;
45 
46  Vector vap = Pos - a, vbp = Pos - b,
47  // vcp = Pos - c,
48  // vdp = Pos - d,
49  vab = b - a, vac = c - a, vad = d - a, vbc = c - b, vbd = d - b;
50 
51  Vector temp_ac_ad(vac[1] * vad[2] - vac[2] * vad[1],
52  vac[2] * vad[0] - vac[0] * vad[2],
53  vac[0] * vad[1] - vac[1] * vad[0]);
54 
55  Vector temp(vbd[1] * vbc[2] - vbd[2] * vbc[1],
56  vbd[2] * vbc[0] - vbd[0] * vbc[2],
57  vbd[0] * vbc[1] - vbd[1] * vbc[0]);
58 
59  Scalar va = (vbp * temp) * 1.f / 6;
60 
61  Scalar vb = (vap * temp_ac_ad) * 1.f / 6;
62 
63  temp = Vector(vad[1] * vab[2] - vad[2] * vab[1],
64  vad[2] * vab[0] - vad[0] * vab[2],
65  vad[0] * vab[1] - vad[1] * vab[0]);
66 
67  Scalar vc = (vap * temp) * 1.f / 6;
68 
69  temp = Vector(vab[1] * vac[2] - vab[2] * vac[1],
70  vab[2] * vac[0] - vab[0] * vac[2],
71  vab[0] * vac[1] - vab[1] * vac[0]);
72 
73  Scalar vd = (vap * temp) * 1.f / 6;
74 
75  Scalar v = (vab * temp_ac_ad) * 1.f / 6;
76 
77  Scalar alpha = va / v, beta = vb / v, gamma = vc / v, delta = vd / v;
79  // std::cout << "barycentric coordinates: (" << alpha << ", " << beta << ", "
80  // << gamma << ", " << delta << ")" << std::endl ;
81 
82  return !((alpha < 0.f) || (beta < 0.f) || (gamma < 0.f) || (delta < 0.f));
83 }
84 
107 template< typename GeometryTraits >
108 static inline
109 typename GeometryTraits::Scalar
110 trilinear_interpolation(const typename GeometryTraits::Point &a,
111  const typename GeometryTraits::Point &b,
112  const typename GeometryTraits::Point &c,
113  const typename GeometryTraits::Point &d,
114  const typename GeometryTraits::Point &Pos,
115  typename GeometryTraits::Scalar vala,
116  typename GeometryTraits::Scalar valb,
117  typename GeometryTraits::Scalar valc,
118  typename GeometryTraits::Scalar vald,
119  const GeometryTraits &gt)
120 {
121  typedef typename GeometryTraits::Vector Vector;
122  typedef typename GeometryTraits::Scalar Scalar;
123 
124  Vector vap = Pos - a, vbp = Pos - b,
125  // vcp = Pos - c,
126  // vdp = Pos - d,
127  vab = b - a, vac = c - a, vad = d - a, vbc = c - b, vbd = d - b;
128 
129  Vector temp_ac_ad(vac[1] * vad[2] - vac[2] * vad[1],
130  vac[2] * vad[0] - vac[0] * vad[2],
131  vac[0] * vad[1] - vac[1] * vad[0]);
132 
133  Vector temp(vbd[1] * vbc[2] - vbd[2] * vbc[1],
134  vbd[2] * vbc[0] - vbd[0] * vbc[2],
135  vbd[0] * vbc[1] - vbd[1] * vbc[0]);
136 
137  Scalar va = (vbp * temp) * 1.f / 6;
138 
139  Scalar vb = (vap * temp_ac_ad) * 1.f / 6;
140 
141  temp = Vector(vad[1] * vab[2] - vad[2] * vab[1],
142  vad[2] * vab[0] - vad[0] * vab[2],
143  vad[0] * vab[1] - vad[1] * vab[0]);
144 
145  Scalar vc = (vap * temp) * 1.f / 6;
146 
147  temp = Vector(vab[1] * vac[2] - vab[2] * vac[1],
148  vab[2] * vac[0] - vab[0] * vac[2],
149  vab[0] * vac[1] - vab[1] * vac[0]);
150 
151  Scalar vd = (vap * temp) * 1.f / 6;
152 
153  Scalar v = (vab * temp_ac_ad) * 1.f / 6;
154 
155  Scalar alpha = va / v, beta = vb / v, gamma = vc / v, delta = vd / v;
157  if(alpha < 0.f)
158  alpha = 0.f;
159  if(alpha > 1.f)
160  alpha = 1.f;
161 
162  if(beta < 0.f)
163  beta = 0.f;
164  if(beta > 1.f)
165  beta = 1.f;
166 
167  if(gamma < 0.f)
168  gamma = 0.f;
169  if(gamma > 1.f)
170  gamma = 1.f;
171 
172  if(delta < 0.f)
173  delta = 0.f;
174  if(delta > 1.f)
175  delta = 1.f;
177  return alpha * vala + beta * valb + gamma * valc + delta * vald;
178 };
179 
202 template< typename GeometryTraits >
203 static inline std::vector< typename GeometryTraits::Scalar >
204 trilinear_interpolation(const typename GeometryTraits::Point &a,
205  const typename GeometryTraits::Point &b,
206  const typename GeometryTraits::Point &c,
207  const typename GeometryTraits::Point &d,
208  const typename GeometryTraits::Point &Pos,
209  const std::vector< typename GeometryTraits::Scalar > &vala,
210  const std::vector< typename GeometryTraits::Scalar > &valb,
211  const std::vector< typename GeometryTraits::Scalar > &valc,
212  const std::vector< typename GeometryTraits::Scalar > &vald,
213  const GeometryTraits &gt)
214 {
215  size_t nbe = vala.size();
216  if(valb.size() < nbe)
217  nbe = valb.size();
218  if(valc.size() < nbe)
219  nbe = valc.size();
220  if(vald.size() < nbe)
221  nbe = vald.size();
223  std::vector< typename GeometryTraits::Scalar > res(nbe, 0.f);
224 
225  for(size_t i = 0; i < nbe; ++i)
226  {
227  res[i] = trilinear_interpolation< GeometryTraits >(
228  a, b, c, d, Pos, vala[i], valb[i], valc[i], vald[i], gt);
229  }
230 
231  return res;
232 };
233 
250 template< typename GeometryTraits >
251 static inline bool
252 is_over_segment(const typename GeometryTraits::Point &a,
253  const typename GeometryTraits::Point &b,
254  const typename GeometryTraits::Point &Pos,
255  typename GeometryTraits::Vector &diff,
256  const GeometryTraits &gt)
257 {
258  typedef typename GeometryTraits::Vector Vector;
259  typedef typename GeometryTraits::Point Point;
260  Vector vab(gt.sub_p(b, a)), vaPos(gt.sub_p(Pos, a)), inter, N, Ntmp;
261  Point interp;
262  inter = gt.cross_product(vab, vaPos);
263  // std::cout << " inter: (" << inter[0] << " " << inter[1] << " " << inter[2]
264  // << ")" << std::endl;
265  Ntmp = gt.cross_product(inter, vab);
266  // std::cout << " Ntmp: (" << Ntmp[0] << " " << Ntmp[1] << " " << Ntmp[2] <<
267  // ")" << std::endl;
268  N = gt.normalize(Ntmp);
269  // std::cout << " N: (" << N[0] << " " << N[1] << " " << N[2] << ")" <<
270  // std::endl;
271  diff = gt.scalar_mult(N,
272  -gt.dot_product(vaPos, N));
273  // "-" to go towards the surface!!
274 
275  interp = Point(Pos[0] + diff[0], Pos[1] + diff[1], Pos[2] + diff[2]);
276  // std::cout << " interp: (" << interp[0] << " " << interp[1] << " " <<
277  // interp[2] << ")" << std::endl;
278  return (gt.dot_product(
279  Vector(interp[0] - a[0], interp[1] - a[1], interp[2] - a[2]),
280  vab) >= 0) &&
281  (gt.dot_product(
282  Vector(b[0] - interp[0], b[1] - interp[1], b[2] - interp[2]),
283  vab) >= 0);
284 }
285 
301 template< typename GeometryTraits >
302 static inline bool
303 is_in_triangle(const typename GeometryTraits::Point &a,
304  const typename GeometryTraits::Point &b,
305  const typename GeometryTraits::Point &c,
306  const typename GeometryTraits::Point &Pos,
307  const GeometryTraits &gt)
308 {
309  typedef typename GeometryTraits::Vector Vector;
310  typedef typename GeometryTraits::Scalar Scalar;
311  // Pos =
312  // we are going to compute the barycentric coordinates of Pos
313  Vector p1p(a[0] - c[0], a[1] - c[1], a[2] - c[2]),
314  p2p(b[0] - c[0], b[1] - c[1], b[2] - c[2]),
315  Posp(Pos[0] - c[0], Pos[1] - c[1], Pos[2] - c[2]);
316 
317  Scalar p1p2 = gt.dot_product(p1p, p1p), p2p2 = gt.dot_product(p2p, p2p),
318  p1pp2p = gt.dot_product(p1p, p2p), p1pPosp = gt.dot_product(p1p, Posp),
319  p2pPosp = gt.dot_product(p2p, Posp);
321  Scalar den = p1p2 * p2p2 - p1pp2p * p1pp2p;
322  Scalar alpha = (p2p2 * p1pPosp - p1pp2p * p2pPosp) / den,
323  beta = (p1p2 * p2pPosp - p1pp2p * p1pPosp) / den, gamma;
324  gamma = 1 - (alpha + beta);
326  return !((alpha < 0.0) || (beta < 0.0) || (gamma < 0.0) || (alpha > 1.0) ||
327  (beta > 1.0) || (gamma > 1.0));
328 };
329 
350 template< typename GeometryTraits >
351 static inline typename GeometryTraits::Scalar
352 bilinear_interpolation(const typename GeometryTraits::Point &a,
353  const typename GeometryTraits::Point &b,
354  const typename GeometryTraits::Point &c,
355  const typename GeometryTraits::Point &Pos,
356  typename GeometryTraits::Scalar vala,
357  typename GeometryTraits::Scalar valb,
358  typename GeometryTraits::Scalar valc,
359  const GeometryTraits &gt)
360 {
361  typedef typename GeometryTraits::Vector Vector;
362  typedef typename GeometryTraits::Scalar Scalar;
363 
364  Vector p1p = a - c, p2p = b - c, Posp = Pos - c;
365  Scalar p1p2 = gt.dot_product(p1p, p1p), p2p2 = gt.dot_product(p2p, p2p),
366  p1pp2p = gt.dot_product(p1p, p2p), p1pPosp = gt.dot_product(p1p, Posp),
367  p2pPosp = gt.dot_product(p2p, Posp);
369  Scalar den = p1p2 * p2p2 - p1pp2p * p1pp2p;
370  Scalar alpha = (p2p2 * p1pPosp - p1pp2p * p2pPosp) / den,
371  beta = (p1p2 * p2pPosp - p1pp2p * p1pPosp) / den, gamma;
372  gamma = 1 - (alpha + beta);
374  if(alpha < 0.f)
375  alpha = 0.f;
376  if(alpha > 1.f)
377  alpha = 1.f;
378 
379  if(beta < 0.f)
380  beta = 0.f;
381  if(beta > 1.f)
382  beta = 1.f;
383 
384  if(gamma < 0.f)
385  gamma = 0.f;
386  if(gamma > 1.f)
387  gamma = 1.f;
389  return alpha * vala + beta * valb + gamma * valc;
390 };
411 template< typename GeometryTraits >
412 static inline std::vector< typename GeometryTraits::Scalar >
413 bilinear_interpolation(
414  const typename GeometryTraits::Point &a,
415  const typename GeometryTraits::Point &b,
416  const typename GeometryTraits::Point &c,
417  const typename GeometryTraits::Point &Pos,
418  const std::vector< typename GeometryTraits::Scalar > &vala,
419  const std::vector< typename GeometryTraits::Scalar > &valb,
420  const std::vector< typename GeometryTraits::Scalar > &valc,
421  const GeometryTraits &gt)
422 {
423  size_t nbe = vala.size();
424  if(valb.size() < nbe)
425  nbe = valb.size();
426  if(valc.size() < nbe)
427  nbe = valc.size();
429  std::vector< typename GeometryTraits::Scalar > res(nbe, 0);
430 
431  for(size_t i = 0; i < nbe; ++i)
432  {
433  res[i] = bilinear_interpolation< GeometryTraits >(
434  a, b, c, Pos, vala[i], valb[i], valc[i], gt);
435  }
436 
437  return res;
438 };
439 } // namespace Geometry
440 } // namespace Operators
441 } // namespace FEVV
Vector
AIFMesh::Vector Vector
Definition: Graph_properties_aif.h:22
Point
AIFMesh::Point Point
Definition: Graph_properties_aif.h:21
FEVV::Operators::Geometry::is_in_tetrahedron
static bool is_in_tetrahedron(const typename GeometryTraits::Point &a, const typename GeometryTraits::Point &b, const typename GeometryTraits::Point &c, const typename GeometryTraits::Point &d, const typename GeometryTraits::Point &Pos, const GeometryTraits &gt)
Predicate fonction to know whether the Pos position is inside the abcd tetrahedron.
Definition: Interpolation.hpp:36
FEVV
Interfaces for plugins These interfaces will be used for different plugins.
Definition: Assert.h:16
FEVV::DataStructures::AIF::AIFVector
Definition: AIFProperties.h:173
FEVV::Operators::Geometry::is_over_segment
static bool is_over_segment(const typename GeometryTraits::Point &a, const typename GeometryTraits::Point &b, const typename GeometryTraits::Point &Pos, typename GeometryTraits::Vector &diff, const GeometryTraits &gt)
Predicate fonction to know whether the Pos position is over the ab segment.
Definition: Interpolation.hpp:252
FEVV::Operators::Geometry::is_in_triangle
static bool is_in_triangle(const typename GeometryTraits::Point &a, const typename GeometryTraits::Point &b, const typename GeometryTraits::Point &c, const typename GeometryTraits::Point &Pos, const GeometryTraits &gt)
Predicate fonction to know whether the projection of point Pos on the triangle plane is in the triang...
Definition: Interpolation.hpp:303
FEVV::Operators::Geometry::trilinear_interpolation
static GeometryTraits::Scalar trilinear_interpolation(const typename GeometryTraits::Point &a, const typename GeometryTraits::Point &b, const typename GeometryTraits::Point &c, const typename GeometryTraits::Point &d, const typename GeometryTraits::Point &Pos, typename GeometryTraits::Scalar vala, typename GeometryTraits::Scalar valb, typename GeometryTraits::Scalar valc, typename GeometryTraits::Scalar vald, const GeometryTraits &gt)
Trilinear interpolation of values at tetrahedron points.
Definition: Interpolation.hpp:110
FEVV::Operators::Geometry::bilinear_interpolation
static GeometryTraits::Scalar bilinear_interpolation(const typename GeometryTraits::Point &a, const typename GeometryTraits::Point &b, const typename GeometryTraits::Point &c, const typename GeometryTraits::Point &Pos, typename GeometryTraits::Scalar vala, typename GeometryTraits::Scalar valb, typename GeometryTraits::Scalar valc, const GeometryTraits &gt)
Bilinear interpolation of values at triangle points.
Definition: Interpolation.hpp:352
FEVV::DataStructures::AIF::AIFPoint
Definition: AIFProperties.h:31
Generated by doxygen 1.8.20