dcv.imgproc.interpolate source code

1 /**
2 Value interpolation module.
3 
4 Copyright: Copyright Relja Ljubobratovic 2016.
5 
6 Authors: Relja Ljubobratovic
7 
8 License: $(LINK3 http://www.boost.org/LICENSE_1_0.txt, Boost Software License - Version 1.0).
9 */
10 module dcv.imgproc.interpolate;
11 
12 import std.traits : isNumeric, isScalarType, isIntegral, allSameType, ReturnType,
13     isFloatingPoint;
14 import std.meta : allSatisfy;
15 
16 import mir.math.common : fastmath;
17 
18 import mir.ndslice.slice;
19 
20 
21 /**
22 Test if given function is proper form for interpolation.
23 */
24 static bool isInterpolationFunc(alias F)()
25 {
26     auto s = [0., 1.].sliced(2);
27     return (__traits(compiles, F(s, 0))); // TODO: check the return type?
28 }
29 
30 /**
31 Test for 1D (vector) interpolation function.
32 */
33 static bool isInterpolationFunc1D(alias F)()
34 {
35     return isInterpolationFunc!F;
36 }
37 
38 /**
39 Test for 2D (matrix) interpolation function.
40 */
41 static bool isInterpolationFunc2D(alias F)()
42 {
43     auto s = [0, 1, 2, 3].sliced(2, 2);
44     return (__traits(compiles, F(s, 3, 3)));
45 }
46 
47 unittest
48 {
49     static assert(isInterpolationFunc!linear);
50     static assert(isInterpolationFunc1D!linear);
51     static assert(isInterpolationFunc2D!linear);
52 }
53 
54 /**
55 Linear interpolation.
56 
57 Params:
58     slice = Input slice which values are interpolated.
59     pos = Position on which slice values are interpolated.
60 
61 Returns:
62     Interpolated resulting value.
63 */
64 pure T linear(SliceKind kind, size_t[] packs, T, P, size_t N)(Slice!(kind, packs, T*) slice, P[N] pos...)
65     if (packs.length == 1)
66 {
67     // TODO: document
68     static assert(N == packs[0], "Interpolation indexing has to be of same dimension as the input slice.");
69 
70     static if (pos.length == 1)
71     {
72         return linearImpl_1(slice, pos[0]);
73     }
74     else static if (pos.length == 2)
75     {
76         return linearImpl_2(slice, pos[0], pos[1]);
77     }
78     else
79     {
80         static assert(0, "Unsupported slice dimension for linear interpolation.");
81     }
82 }
83 
84 unittest
85 {
86     auto arr1 = [0., 1.].sliced(2);
87     assert(linear(arr1, 0.) == 0.);
88     assert(linear(arr1, 1.) == 1.);
89     assert(linear(arr1, 0.1) == 0.1);
90     assert(linear(arr1, 0.5) == 0.5);
91     assert(linear(arr1, 0.9) == 0.9);
92 
93     auto arr1_integral = [0, 10].sliced(2);
94     assert(linear(arr1_integral, 0.) == 0);
95     assert(linear(arr1_integral, 1.) == 10);
96     assert(linear(arr1_integral, 0.1) == 1);
97     assert(linear(arr1_integral, 0.5) == 5);
98     assert(linear(arr1_integral, 0.9) == 9);
99 
100     auto arr2 = [0., 0., 0., 1.].sliced(2, 2);
101     assert(arr2.linear(0.5, 0.5) == 0.25);
102     assert(arr2.linear(0., 0.) == 0.);
103     assert(arr2.linear(1., 1.) == 1.);
104     assert(arr2.linear(1., 0.) == 0.);
105 }
106 
107 private:
108 
109 pure @fastmath auto linearImpl_1(T)(Slice!(SliceKind.contiguous, [1], T*) range, double pos)
110 {
111     import mir.math.common;
112 
113     assert(pos < range.length);
114 
115     if (pos == range.length - 1)
116     {
117         return range[$ - 1];
118     }
119 
120     size_t round = cast(size_t)pos.floor;
121     double weight = pos - cast(double)round;
122 
123     static if (isIntegral!T)
124     {
125         // TODO: is this branch really necessary?
126         auto v1 = cast(double)range[round];
127         auto v2 = cast(double)range[round + 1];
128     }
129     else
130     {
131         auto v1 = range[round];
132         auto v2 = range[round + 1];
133     }
134     return cast(T)(v1 * (1. - weight) + v2 * (weight));
135 }
136 
137 pure @fastmath auto linearImpl_2(SliceKind kind, size_t[] packs, T)(Slice!(kind, packs, T*) range, double pos_x, double pos_y)
138 {
139     import mir.math.common : floor;
140 
141     assert(pos_x < range.length!0 && pos_y < range.length!1);
142 
143     size_t rx = cast(size_t)pos_x.floor;
144     size_t ry = cast(size_t)pos_y.floor;
145     double wx = pos_x - cast(double)rx;
146     double wy = pos_y - cast(double)ry;
147 
148     auto w00 = (1. - wx) * (1. - wy);
149     auto w01 = (wx) * (1. - wy);
150     auto w10 = (1. - wx) * (wy);
151     auto w11 = (wx) * (wy);
152 
153     auto x_end = rx == range.length!0 - 1;
154     auto y_end = ry == range.length!1 - 1;
155 
156     static if (isIntegral!T)
157     {
158         // TODO: (same as in 1D vesion) is this branch really necessary?
159         double v1, v2, v3, v4;
160         v1 = cast(double)range[rx, ry];
161         v2 = cast(double)range[x_end ? rx : rx + 1, ry];
162         v3 = cast(double)range[rx, y_end ? ry : ry + 1];
163         v4 = cast(double)range[x_end ? rx : rx + 1, y_end ? ry : ry + 1];
164     }
165     else
166     {
167         T v1, v2, v3, v4;
168         v1 = range[rx, ry];
169         v2 = range[x_end ? rx : rx + 1, ry];
170         v3 = range[rx, y_end ? ry : ry + 1];
171         v4 = range[x_end ? rx : rx + 1, y_end ? ry : ry + 1];
172     }
173     return cast(T)(v1 * w00 + v2 * w01 + v3 * w10 + v4 * w11);
174 }
175