1 /**
2 inmath.lerpolate
3
4 Authors: David Herberth, Inochi2D Project
5 License: MIT
6 */
7 module inmath.interpolate;
8
9 private {
10 import inmath.linalg : Vector, dot, vec2, vec3, vec4, quat;
11 import inmath.util : isVector, isQuaternion;
12 import inmath.math : almostEqual, acos, sin, sqrt, clamp, PI;
13 import std.conv : to;
14 }
15
16 @safe pure nothrow:
17
18 /// lerpolates linear between two points, also known as lerp.
19 T lerp(T)(T a, T b, float t) {
20 return a * (1 - t) + b * t;
21 }
22
23 /// lerpolates spherical between to vectors or quaternions, also known as slerp.
24 T slerp(T)(T a, T b, float t) if(isVector!T || isQuaternion!T) {
25 static if(isVector!T) {
26 real theta = acos(dot(a, b));
27 } else {
28 real theta = acos(
29 // this is a workaround, acos returning -nan on certain values near +/-1
30 clamp(a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z, -1, 1)
31 );
32 }
33
34 if(almostEqual(theta, 0)) {
35 return a;
36 } else if(almostEqual(theta, PI)) { // 180°?
37 return lerp(a, b, t);
38 } else { // slerp
39 real sintheta = sin(theta);
40 return (sin((1.0-t)*theta)/sintheta)*a + (sin(t*theta)/sintheta)*b;
41 }
42 }
43
44
45 /// Normalized quaternion linear lerpolation.
46 quat nlerp(quat a, quat b, float t) {
47 // TODO: tests
48 float dot = a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z;
49
50 quat result;
51 if(dot < 0) { // Determine the "shortest route"...
52 result = a - (b + a) * t; // use -b instead of b
53 } else {
54 result = a + (b - a) * t;
55 }
56 result.normalize();
57
58 return result;
59 }
60
61 @("vector lerp, slerp")
62 unittest {
63 vec2 v2_1 = vec2(1.0f);
64 vec2 v2_2 = vec2(0.0f);
65 vec3 v3_1 = vec3(1.0f);
66 vec3 v3_2 = vec3(0.0f);
67 vec4 v4_1 = vec4(1.0f);
68 vec4 v4_2 = vec4(0.0f);
69
70 assert(lerp(v2_1, v2_2, 0.5f).vector == [0.5f, 0.5f]);
71 assert(lerp(v2_1, v2_2, 0.0f) == v2_1);
72 assert(lerp(v2_1, v2_2, 1.0f) == v2_2);
73 assert(lerp(v3_1, v3_2, 0.5f).vector == [0.5f, 0.5f, 0.5f]);
74 assert(lerp(v3_1, v3_2, 0.0f) == v3_1);
75 assert(lerp(v3_1, v3_2, 1.0f) == v3_2);
76 assert(lerp(v4_1, v4_2, 0.5f).vector == [0.5f, 0.5f, 0.5f, 0.5f]);
77 assert(lerp(v4_1, v4_2, 0.0f) == v4_1);
78 assert(lerp(v4_1, v4_2, 1.0f) == v4_2);
79
80 assert(slerp(v2_1, v2_2, 0.0).vector == v2_1.vector);
81 assert(slerp(v2_1, v2_2, 1.0).vector == v2_2.vector);
82 assert(slerp(v3_1, v3_2, 0.0).vector == v3_1.vector);
83 assert(slerp(v3_1, v3_2, 1.0).vector == v3_2.vector);
84 assert(slerp(v4_1, v4_2, 0.0).vector == v4_1.vector);
85 assert(slerp(v4_1, v4_2, 1.0).vector == v4_2.vector);
86 }
87
88 @("real lerp")
89 unittest {
90 real r1 = 0.0;
91 real r2 = 1.0;
92 assert(lerp(r1, r2, 0.5f) == 0.5);
93 assert(lerp(r1, r2, 0.0f) == r1);
94 assert(lerp(r1, r2, 1.0f) == r2);
95
96 assert(lerp(0.0, 1.0, 0.5f) == 0.5);
97 assert(lerp(0.0, 1.0, 0.0f) == 0.0);
98 assert(lerp(0.0, 1.0, 1.0f) == 1.0);
99
100 assert(lerp(0.0f, 1.0f, 0.5f) == 0.5f);
101 assert(lerp(0.0f, 1.0f, 0.0f) == 0.0f);
102 assert(lerp(0.0f, 1.0f, 1.0f) == 1.0f);
103 }
104
105 @("quat lerp")
106 unittest {
107 quat q1 = quat(1.0f, 1.0f, 1.0f, 1.0f);
108 quat q2 = quat(0.0f, 0.0f, 0.0f, 0.0f);
109
110 assert(lerp(q1, q2, 0.0f).quaternion == q1.quaternion);
111 assert(lerp(q1, q2, 0.5f).quaternion == [0.5f, 0.5f, 0.5f, 0.5f]);
112 assert(lerp(q1, q2, 1.0f).quaternion == q2.quaternion);
113
114 assert(slerp(q1, q2, 0.0f).quaternion == q1.quaternion);
115 assert(slerp(q1, q2, 1.0f).quaternion == q2.quaternion);
116 }
117
118 /// Nearest lerpolation of two points.
119 T nearest(T)(T x, T y, float t) {
120 if(t < 0.5f) { return x; }
121 else { return y; }
122 }
123
124 @("nearest lerp")
125 unittest {
126 assert(nearest(0.0, 1.0, 0.5f) == 1.0);
127 assert(nearest(0.0, 1.0, 0.4f) == 0.0);
128 assert(nearest(0.0, 1.0, 0.6f) == 1.0);
129 }
130
131 /// Catmull-rom lerpolation between four points.
132 T catmullrom(T)(T p0, T p1, T p2, T p3, float t) {
133 return 0.5f * ((2 * p1) +
134 (-p0 + p2) * t +
135 (2 * p0 - 5 * p1 + 4 * p2 - p3) * t^^2 +
136 (-p0 + 3 * p1 - 3 * p2 + p3) * t^^3);
137 }
138
139 /// Catmull-derivatives of the lerpolation between four points.
140 T dcatmullrom(T)(T p0, T p1, T p2, T p3, float t) {
141 return 0.5f * ((2 * p1) +
142 (-p0 + p2) +
143 2 * (2 * p0 - 5 * p1 + 4 * p2 - p3) * t +
144 3 * (-p0 + 3 * p1 - 3 * p2 + p3) * t^^2);
145 }
146
147 /// Hermite lerpolation (cubic hermite spline).
148 T hermite(T)(T x, T tx, T y, T ty, float t) {
149 float h1 = 2 * t^^3 - 3 * t^^2 + 1;
150 float h2 = -2* t^^3 + 3 * t^^2;
151 float h3 = t^^3 - 2 * t^^2 + t;
152 float h4 = t^^3 - t^^2;
153 return h1 * x + h3 * tx + h2 * y + h4 * ty;
154 }