Hermes
Loading...
Searching...
No Matches
interpolation.h
Go to the documentation of this file.
1
31
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
40namespace hermes::interpolation {
41
46HERMES_DEVICE_CALLABLE static inline f32 smooth(f32 a, f32 b) {
47 return fmaxf(0.f, 1.f - a / (b * b));
48}
53HERMES_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}
61HERMES_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}
71template<typename T>
72HERMES_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}
83HERMES_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}
92template<typename T>
93HERMES_DEVICE_CALLABLE static inline T mix(T x, T y, f32 a) {
94 return x * (1.f - a) + y * a;
95}
102template<typename T>
103HERMES_DEVICE_CALLABLE static T lerp(T t, T a, T b) { return (1. - t) * a + t * b; }
112template<typename T>
113HERMES_DEVICE_CALLABLE static
114inline 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}
133template<typename T>
134HERMES_DEVICE_CALLABLE static
135inline 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}
146template<typename T>
147HERMES_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}
166template<typename T> HERMES_DEVICE_CALLABLE
167T 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}
213template<typename T> HERMES_DEVICE_CALLABLE static
214float 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}
230template<typename T> HERMES_DEVICE_CALLABLE
231static 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}
252template<typename S, typename T>
253HERMES_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::Index2
Holds 2-dimensional integer index coordinates.
Definition index.h:50
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