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