CoCalc -- optical_flow

GitHub Repository: Tetragramm/opencv
Path: blob/master/modules/video/src/opencl/optical_flow_farneback.cl
¹⁶³³⁷ views
1
/*M///////////////////////////////////////////////////////////////////////////////////////
2
//
3
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4
//
5
//  By downloading, copying, installing or using the software you agree to this license.
6
//  If you do not agree to this license, do not download, install,
7
//  copy or use the software.
8
//
9
//
10
//                           License Agreement
11
//                For Open Source Computer Vision Library
12
//
13
// Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved.
14
// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
15
// Third party copyrights are property of their respective owners.
16
//
17
// @Authors
18
//    Sen Liu, swjtuls1987@126.com
19
//
20
// Redistribution and use in source and binary forms, with or without modification,
21
// are permitted provided that the following conditions are met:
22
//
23
//   * Redistribution's of source code must retain the above copyright notice,
24
//     this list of conditions and the following disclaimer.
25
//
26
//   * Redistribution's in binary form must reproduce the above copyright notice,
27
//     this list of conditions and the following disclaimer in the documentation
28
//     and/or other materials provided with the distribution.
29
//
30
//   * The name of the copyright holders may not be used to endorse or promote products
31
//     derived from this software without specific prior written permission.
32
//
33
// This software is provided by the copyright holders and contributors as is and
34
// any express or implied warranties, including, but not limited to, the implied
35
// warranties of merchantability and fitness for a particular purpose are disclaimed.
36
// In no event shall the Intel Corporation or contributors be liable for any direct,
37
// indirect, incidental, special, exemplary, or consequential damages
38
// (including, but not limited to, procurement of substitute goods or services;
39
// loss of use, data, or profits; or business interruption) however caused
40
// and on any theory of liability, whether in contract, strict liability,
41
// or tort (including negligence or otherwise) arising in any way out of
42
// the use of this software, even if advised of the possibility of such damage.
43
//
44
//M*/
45

46

47
#define tx  (int)get_local_id(0)
48
#define ty  get_local_id(1)
49
#define bx  get_group_id(0)
50
#define bdx (int)get_local_size(0)
51

52
#define BORDER_SIZE 5
53
#define MAX_KSIZE_HALF 100
54

55
#ifndef polyN
56
#define polyN 5
57
#endif
58

59
#if USE_DOUBLE
60
#ifdef cl_amd_fp64
61
#pragma OPENCL EXTENSION cl_amd_fp64:enable
62
#elif defined (cl_khr_fp64)
63
#pragma OPENCL EXTENSION cl_khr_fp64:enable
64
#endif
65
#define TYPE double
66
#define VECTYPE double4
67
#else
68
#define TYPE float
69
#define VECTYPE float4
70
#endif
71

72
__kernel void polynomialExpansion(__global __const float * src, int srcStep,
73
                                  __global float * dst, int dstStep,
74
                                  const int rows, const  int cols,
75
                                  __global __const float * c_g,
76
                                  __global __const float * c_xg,
77
                                  __global __const float * c_xxg,
78
                                  __local float * smem,
79
                                  const VECTYPE ig)
80
{
81
    const int y = get_global_id(1);
82
    const int x = bx * (bdx - 2*polyN) + tx - polyN;
83

84
    int xWarped;
85
    __local float *row = smem + tx;
86

87
    if (y < rows && y >= 0)
88
    {
89
        xWarped = min(max(x, 0), cols - 1);
90

91
        row[0] = src[mad24(y, srcStep, xWarped)] * c_g[0];
92
        row[bdx] = 0.f;
93
        row[2*bdx] = 0.f;
94

95
#pragma unroll
96
        for (int k = 1; k <= polyN; ++k)
97
        {
98
            float t0 = src[mad24(max(y - k, 0), srcStep, xWarped)];
99
            float t1 = src[mad24(min(y + k, rows - 1), srcStep, xWarped)];
100

101
            row[0] += c_g[k] * (t0 + t1);
102
            row[bdx] += c_xg[k] * (t1 - t0);
103
            row[2*bdx] += c_xxg[k] * (t0 + t1);
104
        }
105
    }
106

107
    barrier(CLK_LOCAL_MEM_FENCE);
108

109
    if (y < rows && y >= 0 && tx >= polyN && tx + polyN < bdx && x < cols)
110
    {
111
        TYPE b1 = c_g[0] * row[0];
112
        TYPE b3 = c_g[0] * row[bdx];
113
        TYPE b5 = c_g[0] * row[2*bdx];
114
        TYPE b2 = 0, b4 = 0, b6 = 0;
115

116
#pragma unroll
117
        for (int k = 1; k <= polyN; ++k)
118
        {
119
            b1 += (row[k] + row[-k]) * c_g[k];
120
            b4 += (row[k] + row[-k]) * c_xxg[k];
121
            b2 += (row[k] - row[-k]) * c_xg[k];
122
            b3 += (row[k + bdx] + row[-k + bdx]) * c_g[k];
123
            b6 += (row[k + bdx] - row[-k + bdx]) * c_xg[k];
124
            b5 += (row[k + 2*bdx] + row[-k + 2*bdx]) * c_g[k];
125
        }
126

127
        dst[mad24(y, dstStep, xWarped)] = (float)(b3*ig.s0);
128
        dst[mad24(rows + y, dstStep, xWarped)] = (float)(b2*ig.s0);
129
        dst[mad24(2*rows + y, dstStep, xWarped)] = (float)(b1*ig.s1 + b5*ig.s2);
130
        dst[mad24(3*rows + y, dstStep, xWarped)] = (float)(b1*ig.s1 + b4*ig.s2);
131
        dst[mad24(4*rows + y, dstStep, xWarped)] = (float)(b6*ig.s3);
132
    }
133
}
134

135
inline int idx_row_low(const int y, const int last_row)
136
{
137
    return abs(y) % (last_row + 1);
138
}
139

140
inline int idx_row_high(const int y, const int last_row)
141
{
142
    return abs(last_row - abs(last_row - y)) % (last_row + 1);
143
}
144

145
inline int idx_col_low(const int x, const int last_col)
146
{
147
    return abs(x) % (last_col + 1);
148
}
149

150
inline int idx_col_high(const int x, const int last_col)
151
{
152
    return abs(last_col - abs(last_col - x)) % (last_col + 1);
153
}
154

155
inline int idx_col(const int x, const int last_col)
156
{
157
    return idx_col_low(idx_col_high(x, last_col), last_col);
158
}
159

160
__kernel void gaussianBlur(__global const float * src, int srcStep,
161
                           __global float * dst, int dstStep, const int rows, const  int cols,
162
                           __global const float * c_gKer, const int ksizeHalf,
163
                           __local float * smem)
164
{
165
    const int y = get_global_id(1);
166
    const int x = get_global_id(0);
167

168
    __local float *row = smem + ty * (bdx + 2*ksizeHalf);
169

170
    if (y < rows)
171
    {
172
        // Vertical pass
173
        for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx)
174
        {
175
            int xExt = (int)(bx * bdx) + i - ksizeHalf;
176
            xExt = idx_col(xExt, cols - 1);
177
            row[i] = src[mad24(y, srcStep, xExt)] * c_gKer[0];
178
            for (int j = 1; j <= ksizeHalf; ++j)
179
                row[i] += (src[mad24(idx_row_low(y - j, rows - 1), srcStep, xExt)]
180
                           + src[mad24(idx_row_high(y + j, rows - 1), srcStep, xExt)]) * c_gKer[j];
181
        }
182
    }
183

184
    barrier(CLK_LOCAL_MEM_FENCE);
185

186
    if (y < rows && y >= 0 && x < cols && x >= 0)
187
    {
188
        // Horizontal pass
189
        row += tx + ksizeHalf;
190
        float res = row[0] * c_gKer[0];
191
        for (int i = 1; i <= ksizeHalf; ++i)
192
            res += (row[-i] + row[i]) * c_gKer[i];
193

194
        dst[mad24(y, dstStep, x)] = res;
195
    }
196
}
197

198
__kernel void gaussianBlur5(__global const float * src, int srcStep,
199
                            __global float * dst, int dstStep,
200
                            const int rows, const  int cols,
201
                            __global const float * c_gKer, const int ksizeHalf,
202
                            __local float * smem)
203
{
204
    const int y = get_global_id(1);
205
    const int x = get_global_id(0);
206

207
    const int smw = bdx + 2*ksizeHalf; // shared memory "cols"
208
    __local volatile float *row = smem + 5 * ty * smw;
209

210
    if (y < rows)
211
    {
212
        // Vertical pass
213
        for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx)
214
        {
215
            int xExt = (int)(bx * bdx) + i - ksizeHalf;
216
            xExt = idx_col(xExt, cols - 1);
217

218
#pragma unroll
219
            for (int k = 0; k < 5; ++k)
220
                row[k*smw + i] = src[mad24(k*rows + y, srcStep, xExt)] * c_gKer[0];
221

222
            for (int j = 1; j <= ksizeHalf; ++j)
223
#pragma unroll
224
                for (int k = 0; k < 5; ++k)
225
                    row[k*smw + i] +=
226
                        (src[mad24(k*rows + idx_row_low(y - j, rows - 1), srcStep, xExt)] +
227
                         src[mad24(k*rows + idx_row_high(y + j, rows - 1), srcStep, xExt)]) * c_gKer[j];
228
        }
229
    }
230

231
    barrier(CLK_LOCAL_MEM_FENCE);
232

233
    if (y < rows && y >= 0 && x < cols && x >= 0)
234
    {
235
        // Horizontal pass
236

237
        row += tx + ksizeHalf;
238
        float res[5];
239

240
#pragma unroll
241
        for (int k = 0; k < 5; ++k)
242
            res[k] = row[k*smw] * c_gKer[0];
243

244
        for (int i = 1; i <= ksizeHalf; ++i)
245
#pragma unroll
246
            for (int k = 0; k < 5; ++k)
247
                res[k] += (row[k*smw - i] + row[k*smw + i]) * c_gKer[i];
248

249
#pragma unroll
250
        for (int k = 0; k < 5; ++k)
251
            dst[mad24(k*rows + y, dstStep, x)] = res[k];
252
    }
253
}
254
__constant float c_border[BORDER_SIZE + 1] = { 0.14f, 0.14f, 0.4472f, 0.4472f, 0.4472f, 1.f };
255

256
__kernel void updateMatrices(__global const float * flowx, int xStep,
257
                             __global const float * flowy, int yStep,
258
                             const int rows, const int cols,
259
                             __global const float * R0, int R0Step,
260
                             __global const float * R1, int R1Step,
261
                             __global float * M, int mStep)
262
{
263
    const int y = get_global_id(1);
264
    const int x = get_global_id(0);
265

266
    if (y < rows && y >= 0 && x < cols && x >= 0)
267
    {
268
        float dx = flowx[mad24(y, xStep, x)];
269
        float dy = flowy[mad24(y, yStep, x)];
270
        float fx = x + dx;
271
        float fy = y + dy;
272

273
        int x1 = convert_int(floor(fx));
274
        int y1 = convert_int(floor(fy));
275
        fx -= x1;
276
        fy -= y1;
277

278
        float r2, r3, r4, r5, r6;
279

280
        if (x1 >= 0 && y1 >= 0 && x1 < cols - 1 && y1 < rows - 1)
281
        {
282
            float a00 = (1.f - fx) * (1.f - fy);
283
            float a01 = fx * (1.f - fy);
284
            float a10 = (1.f - fx) * fy;
285
            float a11 = fx * fy;
286

287
            r2 = a00 * R1[mad24(y1, R1Step, x1)] +
288
                 a01 * R1[mad24(y1, R1Step, x1 + 1)] +
289
                 a10 * R1[mad24(y1 + 1, R1Step, x1)] +
290
                 a11 * R1[mad24(y1 + 1, R1Step, x1 + 1)];
291

292
            r3 = a00 * R1[mad24(rows + y1, R1Step, x1)] +
293
                 a01 * R1[mad24(rows + y1, R1Step, x1 + 1)] +
294
                 a10 * R1[mad24(rows + y1 + 1, R1Step, x1)] +
295
                 a11 * R1[mad24(rows + y1 + 1, R1Step, x1 + 1)];
296

297
            r4 = a00 * R1[mad24(2*rows + y1, R1Step, x1)] +
298
                 a01 * R1[mad24(2*rows + y1, R1Step, x1 + 1)] +
299
                 a10 * R1[mad24(2*rows + y1 + 1, R1Step, x1)] +
300
                 a11 * R1[mad24(2*rows + y1 + 1, R1Step, x1 + 1)];
301

302
            r5 = a00 * R1[mad24(3*rows + y1, R1Step, x1)] +
303
                 a01 * R1[mad24(3*rows + y1, R1Step, x1 + 1)] +
304
                 a10 * R1[mad24(3*rows + y1 + 1, R1Step, x1)] +
305
                 a11 * R1[mad24(3*rows + y1 + 1, R1Step, x1 + 1)];
306

307
            r6 = a00 * R1[mad24(4*rows + y1, R1Step, x1)] +
308
                 a01 * R1[mad24(4*rows + y1, R1Step, x1 + 1)] +
309
                 a10 * R1[mad24(4*rows + y1 + 1, R1Step, x1)] +
310
                 a11 * R1[mad24(4*rows + y1 + 1, R1Step, x1 + 1)];
311

312
            r4 = (R0[mad24(2*rows + y, R0Step, x)] + r4) * 0.5f;
313
            r5 = (R0[mad24(3*rows + y, R0Step, x)] + r5) * 0.5f;
314
            r6 = (R0[mad24(4*rows + y, R0Step, x)] + r6) * 0.25f;
315
        }
316
        else
317
        {
318
            r2 = r3 = 0.f;
319
            r4 = R0[mad24(2*rows + y, R0Step, x)];
320
            r5 = R0[mad24(3*rows + y, R0Step, x)];
321
            r6 = R0[mad24(4*rows + y, R0Step, x)] * 0.5f;
322
        }
323

324
        r2 = (R0[mad24(y, R0Step, x)] - r2) * 0.5f;
325
        r3 = (R0[mad24(rows + y, R0Step, x)] - r3) * 0.5f;
326

327
        r2 += r4*dy + r6*dx;
328
        r3 += r6*dy + r5*dx;
329

330
        float scale =
331
            c_border[min(x, BORDER_SIZE)] *
332
            c_border[min(y, BORDER_SIZE)] *
333
            c_border[min(cols - x - 1, BORDER_SIZE)] *
334
            c_border[min(rows - y - 1, BORDER_SIZE)];
335

336
        r2 *= scale;
337
        r3 *= scale;
338
        r4 *= scale;
339
        r5 *= scale;
340
        r6 *= scale;
341

342
        M[mad24(y, mStep, x)] = r4*r4 + r6*r6;
343
        M[mad24(rows + y, mStep, x)] = (r4 + r5)*r6;
344
        M[mad24(2*rows + y, mStep, x)] = r5*r5 + r6*r6;
345
        M[mad24(3*rows + y, mStep, x)] = r4*r2 + r6*r3;
346
        M[mad24(4*rows + y, mStep, x)] = r6*r2 + r5*r3;
347
    }
348
}
349

350
__kernel void boxFilter5(__global const float * src, int srcStep,
351
                         __global float * dst, int dstStep,
352
                         const int rows, const  int cols,
353
                         const int ksizeHalf,
354
                         __local float * smem)
355
{
356
    const int y = get_global_id(1);
357
    const int x = get_global_id(0);
358

359
    const float boxAreaInv = 1.f / ((1 + 2*ksizeHalf) * (1 + 2*ksizeHalf));
360
    const int smw = bdx + 2*ksizeHalf; // shared memory "width"
361
    __local float *row = smem + 5 * ty * smw;
362

363
    if (y < rows)
364
    {
365
        // Vertical pass
366
        for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx)
367
        {
368
            int xExt = (int)(bx * bdx) + i - ksizeHalf;
369
            xExt = min(max(xExt, 0), cols - 1);
370

371
#pragma unroll
372
            for (int k = 0; k < 5; ++k)
373
                row[k*smw + i] = src[mad24(k*rows + y, srcStep, xExt)];
374

375
            for (int j = 1; j <= ksizeHalf; ++j)
376
#pragma unroll
377
                for (int k = 0; k < 5; ++k)
378
                    row[k*smw + i] +=
379
                        src[mad24(k*rows + max(y - j, 0), srcStep, xExt)] +
380
                        src[mad24(k*rows + min(y + j, rows - 1), srcStep, xExt)];
381
        }
382
    }
383

384
    barrier(CLK_LOCAL_MEM_FENCE);
385

386
    if (y < rows && y >= 0 && x < cols && x >= 0)
387
    {
388
        // Horizontal pass
389

390
        row += tx + ksizeHalf;
391
        float res[5];
392

393
#pragma unroll
394
        for (int k = 0; k < 5; ++k)
395
            res[k] = row[k*smw];
396

397
        for (int i = 1; i <= ksizeHalf; ++i)
398
#pragma unroll
399
            for (int k = 0; k < 5; ++k)
400
                res[k] += row[k*smw - i] + row[k*smw + i];
401

402
#pragma unroll
403
        for (int k = 0; k < 5; ++k)
404
            dst[mad24(k*rows + y, dstStep, x)] = res[k] * boxAreaInv;
405
    }
406
}
407

408
__kernel void updateFlow(__global const float * M, int mStep,
409
                         __global float * flowx, int xStep,
410
                         __global float * flowy, int yStep,
411
                         const int rows, const int cols)
412
{
413
    const int y = get_global_id(1);
414
    const int x = get_global_id(0);
415

416
    if (y < rows && y >= 0 && x < cols && x >= 0)
417
    {
418
        float g11 = M[mad24(y, mStep, x)];
419
        float g12 = M[mad24(rows + y, mStep, x)];
420
        float g22 = M[mad24(2*rows + y, mStep, x)];
421
        float h1 =  M[mad24(3*rows + y, mStep, x)];
422
        float h2 =  M[mad24(4*rows + y, mStep, x)];
423

424
        float detInv = 1.f / (g11*g22 - g12*g12 + 1e-3f);
425

426
        flowx[mad24(y, xStep, x)] = (g11*h2 - g12*h1) * detInv;
427
        flowy[mad24(y, yStep, x)] = (g22*h1 - g12*h2) * detInv;
428
    }
429
}
430

431
Product

Resources

Company
