dcv.tracking.opticalflow.pyramidflow source code

1 /**
2 Module contains optical flow pyramid implementation.
3 
4 Pyramidal optical flow evaluation in $(BIG DCV) is designed to be a wrapper to actual optical flow algorithm. $(LINK2 #SparesePyramidFlow, Sparse) and 
5 $(LINK2 #DensePyramidFlow, dense) optical flow algorithms have a corresponding utility class which will evaluate the algorithm 
6 in $(LINK3 https://en.wikipedia.org/wiki/Pyramid_(image_processing)#Gaussian_pyramid, pyramid), coarse-to-fine fashion.
7 
8 ----
9 // Evaluate Horn-Schunck method in pyramid.
10 HornSchunckFlow hsFlow = new HornSchunckFlow(props);
11 DensePyramidFlow densePyramid = new DensePyramidFlow(hsFlow, pyramidLevels); 
12 
13 auto flow = densePyramid.evaluate(current, next);
14 ----
15 
16 Copyright: Copyright Relja Ljubobratovic 2016.
17 
18 Authors: Relja Ljubobratovic
19 
20 License: $(LINK3 http://www.boost.org/LICENSE_1_0.txt, Boost Software License - Version 1.0).
21 */
22 module dcv.tracking.opticalflow.pyramidflow;
23 
24 import mir.ndslice.algorithm: each;
25 
26 import dcv.core.utils : emptySlice;
27 import dcv.core.image;
28 import dcv.imgproc.imgmanip : warp, resize;
29 import dcv.tracking.opticalflow.base;
30 
31 import mir.ndslice.allocation: slice;
32 import mir.ndslice.topology: as, flattened, iota;
33 
34 /**
35 Sparse pyramidal optical flow utility class.
36 */
37 class SparsePyramidFlow : SparseOpticalFlow
38 {
39 
40     private SparseOpticalFlow flowAlgorithm;
41     private uint levelCount;
42 
43     this(SparseOpticalFlow flow, uint levels)
44     in
45     {
46         assert(flow !is null);
47         assert(levels > 0);
48     }
49     body
50     {
51         flowAlgorithm = flow;
52         levelCount = levels;
53     }
54 
55     override float[2][] evaluate(inout Image f1, inout Image f2, in float[2][] points,
56             in float[2][] searchRegions, float[2][] flow = null, bool usePrevious = false)
57     in
58     {
59         assert(!f1.empty && !f2.empty && f1.size == f2.size && f1.channels == 1 && f1.depth == f2.depth);
60         assert(points.length == searchRegions.length);
61         if (usePrevious)
62         {
63             assert(flow !is null);
64             assert(points.length == flow.length);
65         }
66     }
67     body
68     {
69         import std.array : uninitializedArray;
70 
71         size_t[2] size = [f1.height, f1.width];
72         const auto pointCount = points.length;
73 
74         // pyramid flow array - each item is double sized flow from the next
75         size_t[2][] flowPyramid;
76         flowPyramid.length = levelCount;
77         flowPyramid[$ - 1] = size.dup;
78 
79         foreach_reverse (i; 0 .. (levelCount - 1))
80         {
81             size[] /= 2;
82             if (size[0] < 1 || size[1] < 1)
83                 throw new Exception("Pyramid downsampling exceeded minimal image size");
84             flowPyramid[i] = size.dup;
85         }
86 
87         auto flowScale = [cast(float)f1.height / cast(float)flowPyramid[0][0],
88             cast(float)f1.width / cast(float)flowPyramid[0][1]];
89 
90         auto lpoints = points.dup;
91         auto lsearchRegions = searchRegions.dup;
92 
93         alias scale = (ref v) { v[0] /= flowScale[0]; v[1] /= flowScale[1]; };
94         lpoints.sliced.each!scale;
95         lsearchRegions.sliced.each!scale;
96 
97         if (usePrevious)
98         {
99             flow.sliced.each!scale;
100         }
101         else
102         {
103             flow = uninitializedArray!(float[2][])(pointCount);
104             flow[] = [0.0f, 0.0f];
105         }
106 
107         auto h = f1.height;
108         auto w = f1.width;
109 
110         Slice!(SliceKind.contiguous, [2], float*) current, next, f1s, f2s;
111         switch (f1.depth) 
112         {
113             case BitDepth.BD_32:
114                 f1s = f1.sliced!float.flattened.sliced(f1.height, f1.width);
115                 f2s = f2.sliced!float.flattened.sliced(f2.height, f2.width);
116                 break;
117             case BitDepth.BD_16:
118                 f1s = f1.sliced!ushort.flattened.sliced(f1.height, f1.width).as!float.slice;
119                 f2s = f2.sliced!ushort.flattened.sliced(f2.height, f2.width).as!float.slice;
120                 break;
121             default:
122                 f1s = f1.sliced!ubyte.flattened.sliced(f1.height, f1.width).as!float.slice;
123                 f2s = f2.sliced!ubyte.flattened.sliced(f2.height, f2.width).as!float.slice;
124         }
125 
126         // calculate pyramid flow
127         foreach (i; 0 .. levelCount)
128         {
129 
130             auto lh = flowPyramid[i][0];
131             auto lw = flowPyramid[i][1];
132 
133             if (lh != h || lw != w)
134             {
135                 current = f1s.resize([lh, lw]);
136                 next = f2s.resize([lh, lw]);
137             }
138             else
139             {
140                 current = f1s;
141                 next = f2s;
142             }
143 
144             flowAlgorithm.evaluate(current.asImage(f1.format), next.asImage(f2.format), lpoints,
145                     lsearchRegions, flow, true);
146 
147             if (i < levelCount - 1)
148             {
149                 alias twice = (ref v) { v[0] += v[0]; v[1] += v[1]; };
150                 flow.sliced.each!twice;
151                 lpoints.sliced.each!twice;
152                 lsearchRegions.sliced.each!twice;
153             }
154         }
155 
156         return flow;
157     }
158 }
159 
160 /**
161 Dense pyramidal optical flow utility class.
162 */
163 class DensePyramidFlow : DenseOpticalFlow
164 {
165 
166     private DenseOpticalFlow flowAlgorithm;
167     private uint levelCount;
168 
169     this(DenseOpticalFlow flow, uint levels)
170     in
171     {
172         assert(flow !is null);
173         assert(levels > 0);
174     }
175     body
176     {
177         flowAlgorithm = flow;
178         levelCount = levels;
179     }
180 
181     override DenseFlow evaluate(inout Image f1, inout Image f2, DenseFlow prealloc = emptySlice!([3],
182             float), bool usePrevious = false)
183     in
184     {
185         assert(prealloc.length!2 == 2);
186         assert(!f1.empty && f1.size == f2.size && f1.depth == f2.depth && f1.depth == BitDepth.BD_8);
187         if (usePrevious)
188         {
189             assert(prealloc.length!0 == f1.height && prealloc.length!1 == f1.width);
190         }
191     }
192     body
193     {
194         size_t[2] size = [f1.height, f1.width];
195         uint level = 0;
196 
197         // pyramid flow array - each item is double sized flow from the next
198         size_t[2][] flowPyramid;
199         flowPyramid.length = levelCount;
200         flowPyramid[$ - 1] = size.dup;
201 
202         DenseFlow flow;
203 
204         foreach_reverse (i; 0 .. (levelCount - 1))
205         {
206             size[] /= 2;
207             if (size[0] < 1 || size[1] < 1)
208                 throw new Exception("Pyramid downsampling exceeded minimal image size");
209             flowPyramid[i] = size.dup;
210         }
211 
212         // allocate flow for each pyramid level
213         if (usePrevious)
214         {
215             flow = prealloc.resize(flowPyramid[0]);
216         }
217         else
218         {
219             flow = new float[flowPyramid[0][0] * flowPyramid[0][1] * 2].sliced(flowPyramid[0][0], flowPyramid[0][1], 2);
220             flow[] = 0.0f;
221         }
222 
223         auto h = f1.height;
224         auto w = f1.width;
225 
226         Slice!(SliceKind.contiguous, [2], float*) current, next, corig, norig;
227         switch (f1.depth) 
228         {
229             case BitDepth.BD_32:
230                 corig = f1.sliced!float.flattened.sliced(f1.height, f1.width);
231                 norig = f2.sliced!float.flattened.sliced(f2.height, f2.width);
232                 break;
233             case BitDepth.BD_16:
234                 corig = f1.sliced!ushort.flattened.sliced(f1.height, f1.width).as!float.slice;
235                 norig = f2.sliced!ushort.flattened.sliced(f2.height, f2.width).as!float.slice;
236                 break;
237             default:
238                 corig = f1.sliced.flattened.sliced(f1.height, f1.width).as!float.slice;
239                 norig = f2.sliced.flattened.sliced(f2.height, f2.width).as!float.slice;
240         }
241 
242         // first flow used as indicator to skip the first warp.
243         bool firstFlow = usePrevious;
244 
245         // calculate pyramid flow
246         foreach (i; 0 .. levelCount)
247         {
248             auto lh = flow.length!0;
249             auto lw = flow.length!1;
250 
251             if (lh != h || lw != w)
252             {
253                 current = corig.resize([lh, lw]);
254                 next = norig.resize([lh, lw]);
255             }
256             else
257             {
258                 current = corig;
259                 next = norig;
260             }
261 
262             if (!firstFlow)
263             {
264                 // warp the image using previous flow, 
265                 // except if this is the first level
266                 // or usePrevious is false.
267                 current = warp(current, flow);
268             }
269 
270             // evaluate the flow algorithm
271             auto lflow = flowAlgorithm.evaluate(current.asImage(f1.format), next.asImage(f2.format));
272 
273             // add flow calculated in this iteration to previous one.
274             flow[] += lflow;
275 
276             if (i < levelCount - 1)
277             {
278                 flow = flow.resize(flowPyramid[i + 1]);
279                 flow[] *= 2.0f;
280             }
281             // assign the first flow indicator to false.
282             firstFlow = false;
283         }
284 
285         return flow;
286     }
287 
288 }
289 
290 // TODO: implement functional tests.
291 version (unittest)
292 {
293 
294     import std.algorithm.iteration : map;
295     import std.array : array;
296     import std.random : uniform;
297 
298     private auto createImage()
299     {
300         return new Image(32, 32, ImageFormat.IF_MONO, BitDepth.BD_8, (32 * 32)
301                 .iota.map!(v => cast(ubyte)uniform(0, 255)).array);
302     }
303 
304     class DummySparseFlow : SparseOpticalFlow
305     {
306         override float[2][] evaluate(inout Image f1, inout Image f2, in float[2][] points,
307                 in float[2][] searchRegions, float[2][] prevflow = null, bool usePrevious = false)
308         {
309             import std.array : uninitializedArray;
310 
311             return uninitializedArray!(float[2][])(points.length);
312         }
313     }
314 
315     class DummyDenseFlow : DenseOpticalFlow
316     {
317         override DenseFlow evaluate(inout Image f1, inout Image f2, DenseFlow prealloc = emptySlice!([3],
318                 float), bool usePrevious = false)
319         {
320             return new float[f1.height * f1.width * 2].sliced(f1.height, f1.width, 2);
321         }
322     }
323 }
324 
325 unittest
326 {
327     SparsePyramidFlow flow = new SparsePyramidFlow(new DummySparseFlow, 3);
328     auto f1 = createImage();
329     auto f2 = createImage();
330     auto p = 10.iota.map!(v => cast(float[2])[cast(float)uniform(0, 2), cast(float)uniform(0, 2)]).array;
331     auto r = 10.iota.map!(v => cast(float[2])[3.0f, 3.0f]).array;
332     auto f = flow.evaluate(f1, f2, p, r);
333     assert(f.length == p.length);
334 }
335 
336 unittest
337 {
338     DensePyramidFlow flow = new DensePyramidFlow(new DummyDenseFlow, 3);
339     auto f1 = createImage();
340     auto f2 = createImage();
341     auto f = flow.evaluate(f1, f2);
342     assert(f.length!0 == f1.height && f.length!1 == f1.width && f.length!2 == 2);
343 }