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