1 /**
2 Module for various utilities used throughout the library.
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.core.utils;
11
12 import std.traits;
13 import std.meta : allSatisfy;
14
15 import mir.ndslice.slice;
16 import mir.ndslice.topology: iota, flattened;
17
18
19 /// Convenience method to return an empty slice - used mainly as default argument in functions in library.
20 static Slice!(SliceKind.contiguous, packs, V*) emptySlice(size_t[] packs, V)() pure @safe nothrow
21 {
22 return Slice!(SliceKind.contiguous, packs, V*)();
23 }
24
25 package(dcv) @nogc pure nothrow
26 {
27 /**
28 Pack and unpack (N, T*) slices to (N-1, T[M]*).
29 */
30 auto staticPack(size_t CH, SliceKind kind, size_t[] packs, T)(Slice!(kind, packs, T*) slice)
31 if (packs.length == 1)
32 {
33 enum N = packs[0];
34 size_t[N-1] shape = slice.shape[0 .. N-1];
35 ptrdiff_t[N-1] strides = [slice.strides[0] / CH, slice.strides[1] / CH];
36 T[CH]* ptr = cast(T[CH]*)slice._iterator;
37 alias Ret = Slice!(kind, [N-1], T[CH]*);
38 return Ret(shape, strides[0 .. Ret.init._strides.length], ptr);
39 }
40
41 /// ditto
42 auto staticUnpack(size_t CH, SliceKind kind, size_t[] packs, T)(Slice!(kind, packs, T[CH]*) slice)
43 if (packs.length == 1)
44 {
45 enum N = packs[0];
46 size_t[N+1] shape = [slice.shape[0], slice.shape[1], CH];
47 ptrdiff_t[N+1] strides = [slice.strides[0] * CH, slice.strides[1] * CH, 1];
48 T* ptr = cast(T*)slice._iterator;
49 alias Ret = Slice!(kind, [N+1], T*);
50 return Ret(shape, strides[0 .. Ret.init._strides.length], ptr);
51 }
52
53 @safe @nogc nothrow pure auto borders(Shape)(Shape shape, size_t ks)
54 in
55 {
56 static assert(Shape.length == 2, "Non-matrix border extraction is not currently supported.");
57 }
58 body
59 {
60 import std.algorithm.comparison : max;
61
62 static struct Range
63 {
64 import mir.ndslice.iterator;
65 Slice!(SliceKind.contiguous, [1], IotaIterator!ptrdiff_t) rows;
66 Slice!(SliceKind.contiguous, [1], IotaIterator!ptrdiff_t) cols;
67 }
68
69 size_t kh = max(size_t(1), ks / 2);
70
71 Range[4] borders = [
72 Range(iota(shape[0]), iota(kh)),
73 Range(iota(shape[0]), iota([kh], shape[1] - kh)),
74 Range(iota(kh), iota(shape[1])),
75 Range(iota([kh], shape[0] - kh), iota(shape[1])),
76 ];
77
78 return borders;
79 }
80 }
81
82 /**
83 Clip value by it's value range.
84
85 Params:
86 v = input value, of the input type
87
88 Returns:
89 Clipped value of the output type.
90 */
91 static pure nothrow @safe T clip(T, V)(V v) if (isNumeric!V && isNumeric!T)
92 {
93 import std.traits : isFloatingPoint;
94
95 static if (is(T == V))
96 {
97 return v;
98 }
99 else
100 {
101 if (v <= T.min)
102 return T.min;
103 else if (v >= T.max)
104 return T.max;
105 return cast(T)v;
106 }
107 }
108
109 pure nothrow @safe unittest
110 {
111 int value = 30;
112 assert(clip!int(value) == cast(int)30);
113 }
114
115 pure nothrow @safe unittest
116 {
117 int max_ubyte_value = cast(int)ubyte.max + 1;
118 int min_ubyte_value = cast(int)ubyte.min - 1;
119 assert(clip!ubyte(max_ubyte_value) == cast(ubyte)255);
120 assert(clip!ubyte(min_ubyte_value) == cast(ubyte)0);
121 }
122
123 import std.compiler;
124
125 static if (__VERSION__ >= 2071)
126 { // due to undefined bug in flattened in dmd 2.070.0
127
128 /**
129 * Merge multiple slices into one.
130 *
131 * By input of multiple Slice!(kind, [N], T*) objects, produces one Slice!(kind, [N+1], T*)
132 * object, where length of last dimension is number of input slices. Values
133 * of input slices' elements are copied to resulting slice, where [..., i] element
134 * of j-th input slice is copied to [..., i, j] element of output slice.
135 *
136 * e.g. If three single channel images (Slice!(kind, [2], Iterator)) are merged, output will be
137 * a three channel image (Slice!(kind, [3], Iterator)).
138 *
139 * Params:
140 * slices = Input slices. All must by Slice object with same input template parameters.
141 *
142 * Returns:
143 * For input of n Slice!(Contiguous, [N], T*) objects, outputs Slice!(kind, [N+1], T*) object, where
144 * last dimension size is n.
145 */
146 pure auto merge(Slices...)(Slices slices)
147 if (Slices.length > 0 && isSlice!(Slices[0]) && allSameType!Slices)
148 {
149 import std.algorithm.iteration : map;
150 import std.array : uninitializedArray;
151
152 alias ElementRange = typeof(slices[0].flattened);
153 alias T = typeof(slices[0].flattened.front);
154
155 immutable D = slices[0].shape.length;
156 const auto length = slices[0].elementsCount;
157
158 auto data = uninitializedArray!(T[])(length * slices.length);
159 ElementRange[slices.length] elRange;
160
161 foreach (i, v; slices)
162 {
163 elRange[i] = v.flattened;
164 }
165
166 auto i = 0;
167 foreach (e; 0 .. length)
168 {
169 foreach (ecol; elRange)
170 {
171 data[i++] = ecol[e];
172 }
173 }
174
175 size_t[D + 1] newShape;
176 newShape[0 .. D] = slices[0].shape[0 .. D];
177 newShape[D] = slices.length;
178
179 return data.sliced(newShape);
180 }
181
182 unittest
183 {
184 auto s1 = [1, 2, 3].sliced(3);
185 auto s2 = [4, 5, 6].sliced(3);
186 auto m = merge(s1, s2);
187 assert(m == [1, 4, 2, 5, 3, 6].sliced(3, 2));
188 }
189
190 unittest
191 {
192 auto s1 = [1, 2, 3, 4].sliced(2, 2);
193 auto s2 = [5, 6, 7, 8].sliced(2, 2);
194 auto m = merge(s1, s2);
195 assert(m == [1, 5, 2, 6, 3, 7, 4, 8].sliced(2, 2, 2));
196 }
197 }
198
199 /// Check if given function can perform boundary condition test.
200 template isBoundaryCondition(alias bc)
201 {
202 enum bool isBoundaryCondition = true;
203 }
204
205 nothrow @nogc pure
206 {
207
208 /// $(LINK2 https://en.wikipedia.org/wiki/Neumann_boundary_condition, Neumann) boundary condition test
209 auto neumann(SliceKind kind, size_t[] packs, Iterator, size_t N)(Slice!(kind, packs, Iterator) tensor, size_t[N] indices...)
210 {
211 static assert(packs[0] == N, "Invalid index dimension");
212 return tensor.bcImpl!_neumann(indices);
213 }
214
215 /// $(LINK2 https://en.wikipedia.org/wiki/Periodic_boundary_conditions,Periodic) boundary condition test
216 auto periodic(SliceKind kind, size_t[] packs, Iterator, size_t N)(Slice!(kind, packs, Iterator) tensor, size_t[N] indices...)
217 {
218 static assert(packs[0] == N, "Invalid index dimension");
219 return tensor.bcImpl!_periodic(indices);
220 }
221
222 /// Symmetric boundary condition test
223 auto symmetric(SliceKind kind, size_t[] packs, Iterator, size_t N)(Slice!(kind, packs, Iterator) tensor, size_t[N] indices...)
224 {
225 static assert(packs[0] == N, "Invalid index dimension");
226 return tensor.bcImpl!_symmetric(indices);
227 }
228
229 }
230
231 /// Alias for generalized boundary condition test function.
232 template BoundaryConditionTest(SliceKind kind, size_t[] packs, T, size_t N)
233 {
234 alias BoundaryConditionTest = ref T function(ref Slice!(kind, packs, T*) slice, size_t[N] indices...);
235 }
236
237 unittest
238 {
239 import mir.ndslice.allocation;
240
241 static assert(isBoundaryCondition!neumann);
242 static assert(isBoundaryCondition!periodic);
243 static assert(isBoundaryCondition!symmetric);
244
245 /*
246 * [0, 1, 2,
247 * 3, 4, 5,
248 * 6, 7, 8]
249 */
250 auto s = iota(3, 3).slice;
251
252 assert(s.neumann(-1, -1) == s[0, 0]);
253 assert(s.neumann(0, -1) == s[0, 0]);
254 assert(s.neumann(-1, 0) == s[0, 0]);
255 assert(s.neumann(0, 0) == s[0, 0]);
256 assert(s.neumann(2, 2) == s[2, 2]);
257 assert(s.neumann(3, 3) == s[2, 2]);
258 assert(s.neumann(1, 3) == s[1, 2]);
259 assert(s.neumann(3, 1) == s[2, 1]);
260
261 assert(s.symmetric(-1, -1) == s[1, 1]);
262 assert(s.symmetric(-2, -2) == s[2, 2]);
263 assert(s.symmetric(0, 0) == s[0, 0]);
264 assert(s.symmetric(2, 2) == s[2, 2]);
265 assert(s.symmetric(3, 3) == s[1, 1]);
266 }
267
268 private:
269
270 nothrow @nogc pure:
271
272 auto bcImpl(alias bcfunc, SliceKind kind, size_t[] packs, Iterator, size_t N)(Slice!(kind, packs, Iterator) tensor, size_t[N] indices...)
273 {
274 static if (N == 1)
275 {
276 return tensor[bcfunc(cast(int)indices[0], cast(int)tensor.length)];
277 }
278 else static if (N == 2)
279 {
280 return tensor[bcfunc(cast(int)indices[0], cast(int)tensor.length!0),
281 bcfunc(cast(int)indices[1], cast(int)tensor.length!1)];
282 }
283 else static if (N == 3)
284 {
285 return tensor[bcfunc(cast(int)indices[0], cast(int)tensor.length!0),
286 bcfunc(cast(int)indices[1], cast(int)tensor.length!1), bcfunc(cast(int)indices[2],
287 cast(int)tensor.length!2)];
288 }
289 else
290 {
291 foreach (i, ref id; indices)
292 {
293 id = bcfunc(cast(int)id, cast(int)tensor._lengths[i]);
294 }
295 return tensor[indices];
296 }
297 }
298
299 int _neumann(int x, int nx)
300 {
301 if (x < 0)
302 {
303 x = 0;
304 }
305 else if (x >= nx)
306 {
307 x = nx - 1;
308 }
309
310 return x;
311 }
312
313 int _periodic(int x, int nx)
314 {
315 if (x < 0)
316 {
317 const int n = 1 - cast(int)(x / (nx + 1));
318 const int ixx = x + n * nx;
319
320 x = ixx % nx;
321 }
322 else if (x >= nx)
323 {
324 x = x % nx;
325 }
326
327 return x;
328 }
329
330 int _symmetric(int x, int nx)
331 {
332 if (x < 0)
333 {
334 const int borde = nx - 1;
335 const int xx = -x;
336 const int n = cast(int)(xx / borde) % 2;
337
338 if (n)
339 x = borde - (xx % borde);
340 else
341 x = xx % borde;
342 }
343 else if (x >= nx)
344 {
345 const int borde = nx - 1;
346 const int n = cast(int)(x / borde) % 2;
347
348 if (n)
349 x = borde - (x % borde);
350 else
351 x = x % borde;
352 }
353 return x;
354 }