Hermes
interpolation.h
Go to the documentation of this file.
1 
32 #ifndef HERMES_HERMES_NUMERIC_INTERPOLATION_H
33 #define HERMES_HERMES_NUMERIC_INTERPOLATION_H
34 
35 #include <algorithm>
36 #include <hermes/common/size.h>
37 #include <hermes/numeric/numeric.h>
38 #include <hermes/geometry/point.h>
39 
40 namespace hermes::interpolation {
41 
46 HERMES_DEVICE_CALLABLE static inline f32 smooth(f32 a, f32 b) {
47  return fmaxf(0.f, 1.f - a / (b * b));
48 }
53 HERMES_DEVICE_CALLABLE static inline f32 sharpen(const f32 &r2, const f32 &h) {
54  return fmaxf(h * h / fmaxf(r2, static_cast<f32>(1.0e-5)) - 1.0f, 0.0f);
55 }
61 HERMES_DEVICE_CALLABLE static f32 smoothStep(f32 a, f32 b, f32 v) {
62  f32 t = Numbers::clamp((v - a) / (b - a), 0.f, 1.f);
63  return t * t * (3.f - 2 * t);
64 }
71 template<typename T>
72 HERMES_DEVICE_CALLABLE static Vector2<T> smoothStep(f32 a, f32 b, const Vector2<T> &v) {
73  Vector2<T> t = {
74  Numbers::clamp((v.x - a) / (b - a), 0.f, 1.f),
75  Numbers::clamp((v.y - a) / (b - a), 0.f, 1.f)};
76  return t * t * (Vector2<T>(3.f, 3.f) - T(2) * t);
77 }
83 HERMES_DEVICE_CALLABLE static inline f32 linearStep(f32 v, f32 a, f32 b) {
84  return Numbers::clamp((v - a) / (b - a), 0.f, 1.f);
85 }
92 template<typename T>
93 HERMES_DEVICE_CALLABLE static inline T mix(T x, T y, f32 a) {
94  return x * (1.f - a) + y * a;
95 }
102 template<typename T>
103 HERMES_DEVICE_CALLABLE static T lerp(T t, T a, T b) { return (1. - t) * a + t * b; }
112 template<typename T>
113 HERMES_DEVICE_CALLABLE static
114 inline T bilerp(T x, T y, const T &f00, const T &f10, const T &f11,
115  const T &f01) {
116  return lerp(y, lerp(x, f00, f10), lerp(x, f01, f11));
117 }
133 template<typename T>
134 HERMES_DEVICE_CALLABLE static
135 inline T trilerp(T tx, T ty, T tz, const T &f000, const T &f100, const T &f010,
136  const T &f110, const T &f001, const T &f101, const T &f011,
137  const T &f111) {
138  return lerp(bilerp(f000, f100, f010, f110, tx, ty),
139  bilerp(f001, f101, f011, f111, tx, ty), tz);
140 }
146 template<typename T>
147 HERMES_DEVICE_CALLABLE static inline T nearest(T t, const T &a, const T &b) {
148  return (t < static_cast<T>(0.5)) ? a : b;
149 }
166 template<typename T> HERMES_DEVICE_CALLABLE
167 T monotonicCubicInterpolate(T fkm1, T fk, T fkp1, T fkp2, T tmtk) {
168  double Dk = fkp1 - fk;
169  double dk = (fkp1 - fkm1) * 0.5;
170  double dkp1 = (fkp2 - fk) * 0.5;
171  if (fabsf(Dk) < 1e-12)
172  dk = dkp1 = 0.0;
173  else {
174  if (Numbers::sign(dk) != Numbers::sign(Dk))
175  // dk = 0;
176  dk *= -1;
177  if (Numbers::sign(dkp1) != Numbers::sign(Dk))
178  // dkp1 = 0;
179  dkp1 *= -1;
180  }
181  double a0 = fk;
182  double a1 = dk;
183  double a2 = 3 * Dk - 2 * dk - dkp1;
184  double a3 = dk + dkp1 - 2 * Dk;
185  T ans = a3 * tmtk * tmtk * tmtk + a2 * tmtk * tmtk + a1 * tmtk + a0;
186 #ifdef HERMES_DEVICE_ENABLED
187  T m = fminf(fkm1, fminf(fk, fminf(fkp1, fkp2)));
188  T M = fmaxf(fkm1, fmaxf(fk, fmaxf(fkp1, fkp2)));
189  return fminf(M, fmaxf(m, ans));
190 #else
191  T m = std::min(fkm1, std::min(fk, std::min(fkp1, fkp2)));
192  T M = std::max(fkm1, std::max(fk, std::max(fkp1, fkp2)));
193  return std::min(M, std::max(m, ans));
194 #endif
195 }
213 template<typename T> HERMES_DEVICE_CALLABLE static
214 float monotonicCubicInterpolate(T f[4][4], const point2 &t) {
215  T v[4];
216  for (int d = 0; d < 4; d++)
217  v[d] = monotonicCubicInterpolate(f[0][d], f[1][d], f[2][d], f[3][d], t.x);
218  return monotonicCubicInterpolate(v[0], v[1], v[2], v[3], t.y);
219 }
230 template<typename T> HERMES_DEVICE_CALLABLE
231 static float monotonicCubicInterpolate(T f[4][4][4], const point3 &t) {
232  float v[4][4];
233  for (int dz = -1, iz = 0; dz <= 2; dz++, iz++)
234  for (int dy = -1, iy = 0; dy <= 2; dy++, iy++)
235  v[iy][iz] = monotonicCubicInterpolate(
236  f[0][dy + 1][dz + 1], f[1][dy + 1][dz + 1], f[2][dy + 1][dz + 1],
237  f[3][dy + 1][dz + 1], t.x);
238  float vv[4];
239  for (int d = 0; d < 4; d++)
240  vv[d] = monotonicCubicInterpolate(v[0][d], v[1][d], v[2][d], v[3][d], t.y);
241  return monotonicCubicInterpolate(vv[0], vv[1], vv[2], vv[3], t.z);
242 }
252 template<typename S, typename T>
253 HERMES_DEVICE_CALLABLE static inline S catmullRomSpline(
254  const S &f0, const S &f1, const S &f2, const S &f3, T f) {
255  S d1 = (f2 - f0) / 2;
256  S d2 = (f3 - f1) / 2;
257  S D1 = f2 - f1;
258 
259  S a3 = d1 + d2 - 2 * D1;
260  S a2 = 3 * D1 - 2 * d1 - d2;
261  S a1 = d1;
262  S a0 = f1;
263 
264  return a3 * CUBE(f) + a2 * SQR(f) + a1 * f + a0;
265 }
266 //template<typename T>
267 //HERMES_DEVICE_CALLABLE static inline T bilinearInterpolation(
268 // T f00, T f10, T f11, T f01, T x, T y) {
269 // return f00 * (1.0 - x) * (1.0 - y) + f10 * x * (1.0 - y) +
270 // f01 * (1.0 - x) * y + f11 * x * y;
271 //}
272 //template<typename T>
273 //HERMES_DEVICE_CALLABLE static inline T cubicInterpolate(T p[4], T x) {
274 // return p[1] + 0.5 * x *
275 // (p[2] - p[0] +
276 // x * (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] +
277 // x * (3.0 * (p[1] - p[2]) + p[3] - p[0])));
278 //}
279 //template<typename T> inline T bicubicInterpolate(T p[4][4], T x, T y) {
280 // T arr[4];
281 // arr[0] = cubicInterpolate(p[0], y);
282 // arr[1] = cubicInterpolate(p[1], y);
283 // arr[2] = cubicInterpolate(p[2], y);
284 // arr[3] = cubicInterpolate(p[3], y);
285 // return cubicInterpolate(arr, x);
286 //}
287 //template<typename T>
288 //inline T trilinearInterpolate(const float *p, T ***data, T b,
289 // const size3 &dimensions) {
290 // i64 dim[3] = {dimensions[0], dimensions[1], dimensions[2]};
291 // i32 i0 = p[0], j0 = p[1], k0 = p[2];
292 // i32 i1 = p[0] + 1, j1 = p[1] + 1, k1 = p[2] + 1;
293 // float x = p[0] - i0;
294 // float y = p[1] - j0;
295 // float z = p[2] - k0;
296 // T v000 = (i0 < 0 || j0 < 0 || k0 < 0 || i0 >= dim[0] ||
297 // j0 >= dim[1] || k0 >= dim[2])
298 // ? b
299 // : data[i0][j0][k0];
300 // T v001 = (i0 < 0 || j0 < 0 || k1 < 0 || i0 >= dim[0] ||
301 // j0 >= dim[1] || k1 >= dim[2])
302 // ? b
303 // : data[i0][j0][k1];
304 // T v010 = (i0 < 0 || j1 < 0 || k0 < 0 || i0 >= dim[0] ||
305 // j1 >= dim[1] || k0 >= dim[2])
306 // ? b
307 // : data[i0][j1][k0];
308 // T v011 = (i0 < 0 || j1 < 0 || k1 < 0 || i0 >= dim[0] ||
309 // j1 >= dim[1] || k1 >= dim[2])
310 // ? b
311 // : data[i0][j1][k1];
312 // T v100 = (i1 < 0 || j0 < 0 || k0 < 0 || i1 >= dim[0] ||
313 // j0 >= dim[1] || k0 >= dim[2])
314 // ? b
315 // : data[i1][j0][k0];
316 // T v101 = (i1 < 0 || j0 < 0 || k1 < 0 || i1 >= dim[0] ||
317 // j0 >= dim[1] || k1 >= dim[2])
318 // ? b
319 // : data[i1][j0][k1];
320 // T v110 = (i1 < 0 || j1 < 0 || k0 < 0 || i1 >= dim[0] ||
321 // j1 >= dim[1] || k0 >= dim[2])
322 // ? b
323 // : data[i1][j1][k0];
324 // T v111 = (i1 < 0 || j1 < 0 || k1 < 0 || i1 >= dim[0] ||
325 // j1 >= dim[1] || k1 >= dim[2])
326 // ? b
327 // : data[i1][j1][k1];
328 // return v000 * (1.f - x) * (1.f - y) * (1.f - z) +
329 // v100 * x * (1.f - y) * (1.f - z) + v010 * (1.f - x) * y * (1.f - z) +
330 // v110 * x * y * (1.f - z) + v001 * (1.f - x) * (1.f - y) * z +
331 // v101 * x * (1.f - y) * z + v011 * (1.f - x) * y * z + v111 * x * y * z;
332 //}
333 //template<typename T> inline T tricubicInterpolate(const float *p, T ***data) {
334 // int x, y, z;
335 // int i, j, k;
336 // float dx, dy, dz;
337 // float u[4], v[4], w[4];
338 // T r[4], q[4];
339 // T vox = T(0);
340 //
341 // x = (int) p[0], y = (int) p[1], z = (int) p[2];
342 // dx = p[0] - (float) x, dy = p[1] - (float) y, dz = p[2] - (float) z;
343 //
344 // u[0] = -0.5f * Numbers::cube(dx) + Numbers::sqr(dx) - 0.5 * dx;
345 // u[1] = 1.5f * Numbers::cube(dx) - 2.5 * Numbers::sqr(dx) + 1;
346 // u[2] = -1.5f * Numbers::cube(dx) + 2 * Numbers::sqr(dx) + 0.5 * dx;
347 // u[3] = 0.5f * Numbers::cube(dx) - 0.5 * Numbers::sqr(dx);
348 //
349 // v[0] = -0.5 * Numbers::cube(dy) + Numbers::sqr(dy) - 0.5 * dy;
350 // v[1] = 1.5 * Numbers::cube(dy) - 2.5 * Numbers::sqr(dy) + 1;
351 // v[2] = -1.5 * Numbers::cube(dy) + 2 * Numbers::sqr(dy) + 0.5 * dy;
352 // v[3] = 0.5 * Numbers::cube(dy) - 0.5 * Numbers::sqr(dy);
353 //
354 // w[0] = -0.5 * Numbers::cube(dz) + Numbers::sqr(dz) - 0.5 * dz;
355 // w[1] = 1.5 * Numbers::cube(dz) - 2.5 * Numbers::sqr(dz) + 1;
356 // w[2] = -1.5 * Numbers::cube(dz) + 2 * Numbers::sqr(dz) + 0.5 * dz;
357 // w[3] = 0.5 * Numbers::cube(dz) - 0.5 * Numbers::sqr(dz);
358 //
359 // int ijk[3] = {x - 1, y - 1, z - 1};
360 // for (k = 0; k < 4; k++) {
361 // q[k] = 0;
362 // for (j = 0; j < 4; j++) {
363 // r[j] = 0;
364 // for (i = 0; i < 4; i++) {
365 // r[j] += u[i] * data[ijk[0]][ijk[1]][ijk[2]];
366 // ijk[0]++;
367 // }
368 // q[k] += v[j] * r[j];
369 // ijk[0] = x - 1;
370 // ijk[1]++;
371 // }
372 // vox += w[k] * q[k];
373 // ijk[0] = x - 1;
374 // ijk[1] = y - 1;
375 // ijk[2]++;
376 // }
377 // return (vox < T(0) ? T(0.0) : vox);
378 //}
379 //template<typename T>
380 //T tricubicInterpolate(const float *p, T ***data, T b, const int dimensions[3]) {
381 // int x, y, z;
382 // int i, j, k;
383 // float dx, dy, dz;
384 // float u[4], v[4], w[4];
385 // T r[4], q[4];
386 // T vox = T(0);
387 //
388 // x = (int) p[0], y = (int) p[1], z = (int) p[2];
389 // dx = p[0] - (float) x, dy = p[1] - (float) y, dz = p[2] - (float) z;
390 //
391 // u[0] = -0.5 * Numbers::cube(dx) + Numbers::sqr(dx) - 0.5 * dx;
392 // u[1] = 1.5 * Numbers::cube(dx) - 2.5 * Numbers::sqr(dx) + 1;
393 // u[2] = -1.5 * Numbers::cube(dx) + 2 * Numbers::sqr(dx) + 0.5 * dx;
394 // u[3] = 0.5 * Numbers::cube(dx) - 0.5 * Numbers::sqr(dx);
395 //
396 // v[0] = -0.5 * Numbers::cube(dy) + Numbers::sqr(dy) - 0.5 * dy;
397 // v[1] = 1.5 * Numbers::cube(dy) - 2.5 * Numbers::sqr(dy) + 1;
398 // v[2] = -1.5 * Numbers::cube(dy) + 2 * Numbers::sqr(dy) + 0.5 * dy;
399 // v[3] = 0.5 * Numbers::cube(dy) - 0.5 * Numbers::sqr(dy);
400 //
401 // w[0] = -0.5 * Numbers::cube(dz) + Numbers::sqr(dz) - 0.5 * dz;
402 // w[1] = 1.5 * Numbers::cube(dz) - 2.5 * Numbers::sqr(dz) + 1;
403 // w[2] = -1.5 * Numbers::cube(dz) + 2 * Numbers::sqr(dz) + 0.5 * dz;
404 // w[3] = 0.5 * Numbers::cube(dz) - 0.5 * Numbers::sqr(dz);
405 //
406 // int ijk[3] = {x - 1, y - 1, z - 1};
407 // for (k = 0; k < 4; k++) {
408 // q[k] = 0;
409 // for (j = 0; j < 4; j++) {
410 // r[j] = 0;
411 // for (i = 0; i < 4; i++) {
412 // if (ijk[0] < 0 || ijk[0] >= dimensions[0] || ijk[1] < 0 ||
413 // ijk[1] >= dimensions[1] || ijk[2] < 0 || ijk[2] >= dimensions[2])
414 // r[j] += u[i] * b;
415 // else
416 // r[j] += u[i] * data[ijk[0]][ijk[1]][ijk[2]];
417 // ijk[0]++;
418 // }
419 // q[k] += v[j] * r[j];
420 // ijk[0] = x - 1;
421 // ijk[1]++;
422 // }
423 // vox += w[k] * q[k];
424 // ijk[0] = x - 1;
425 // ijk[1] = y - 1;
426 // ijk[2]++;
427 // }
428 // return (vox < T(0) ? T(0.0) : vox);
429 //}
430 
431 } // namespace hermes
432 
433 #endif //HERMES_HERMES_NUMERIC_INTERPOLATION_H
434 
hermes::Point2::y
T y
1-th component
Definition: point.h:140
hermes::Point2::x
T x
0-th component
Definition: point.h:139
HERMES_DEVICE_CALLABLE
#define HERMES_DEVICE_CALLABLE
Specifies that the function can be called from both host and device sides.
Definition: defs.h:45
f32
float f32
32 bit size floating point type
Definition: defs.h:78
hermes::interpolation::monotonicCubicInterpolate
HERMES_DEVICE_CALLABLE T monotonicCubicInterpolate(T fkm1, T fk, T fkp1, T fkp2, T tmtk)
Computes 1-dimension monotonic cubic interpolation.
Definition: interpolation.h:167
numeric.h
Number functions.
point.h
Geometric point classes.
size.h
Set of multi-dimensional size representations.
hermes::Numbers::clamp
static HERMES_DEVICE_CALLABLE T clamp(const T &n, const T &l, const T &u)
Clamps value to closed interval.
Definition: numeric.h:176
hermes::Numbers::sign
static HERMES_DEVICE_CALLABLE int sign(T a)
Computes sign.
Definition: numeric.h:202