Coverage Report

Created: 2021-04-14 19:30

/home/liu/buildslave/linux-x64-runtests/build/lib/ccv_classic.c
Line
Count
Source (jump to first uncovered line)
1
#include "ccv.h"
2
#include "ccv_internal.h"
3
4
void ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int b_type, int sbin, int size)
5
0
{
6
0
  assert(a->rows >= size && a->cols >= size && (4 + sbin * 3) <= CCV_MAX_CHANNEL);
7
0
  int rows = a->rows / size;
8
0
  int cols = a->cols / size;
9
0
  b_type = (CCV_GET_DATA_TYPE(b_type) == CCV_64F) ? CCV_64F | (4 + sbin * 3) : CCV_32F | (4 + sbin * 3);
10
0
  ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_hog(%d,%d)", sbin, size), a->sig, CCV_EOF_SIGN);
11
0
  ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, rows, cols, CCV_64F | CCV_32F | (4 + sbin * 3), b_type, sig);
12
0
  ccv_object_return_if_cached(, db);
13
0
  ccv_dense_matrix_t* ag = 0;
14
0
  ccv_dense_matrix_t* mg = 0;
15
0
  ccv_gradient(a, &ag, 0, &mg, 0, 1, 1);
16
0
  float* agp = ag->data.f32;
17
0
  float* mgp = mg->data.f32;
18
0
  int i, j, k, ch = CCV_GET_CHANNEL(a->type);
19
0
  ccv_dense_matrix_t* cn = ccv_dense_matrix_new(rows, cols, CCV_GET_DATA_TYPE(db->type) | (sbin * 2), 0, 0);
20
0
  ccv_dense_matrix_t* ca = ccv_dense_matrix_new(rows, cols, CCV_GET_DATA_TYPE(db->type) | CCV_C1, 0, 0);
21
0
  ccv_zero(cn);
22
0
  // normalize sbin direction-sensitive and sbin * 2 insensitive over 4 normalization factor
23
0
  // accumulating them over sbin * 2 + sbin + 4 channels
24
0
  // TNA - truncation - normalization - accumulation
25
0
#define TNA(_for_type, idx, a, b, c, d) \
26
0
  { \
27
0
    _for_type norm = 1.0 / sqrt(cap[a] + cap[b] + cap[c] + cap[d] + 1e-4); \
28
0
    for (k = 0; k < sbin * 2; k++) \
29
0
    { \
30
0
      _for_type v = 0.5 * ccv_min(cnp[k] * norm, 0.2); \
31
0
      dbp[4 + sbin + k] += v; \
32
0
      dbp[idx] += v; \
33
0
    } \
34
0
    dbp[idx] *= 0.2357; \
35
0
    for (k = 0; k < sbin; k++) \
36
0
    { \
37
0
      _for_type v = 0.5 * ccv_min((cnp[k] + cnp[k + sbin]) * norm, 0.2); \
38
0
      dbp[4 + k] += v; \
39
0
    } \
40
0
  }
41
0
#define for_block(_, _for_type) \
42
0
  _for_type* cnp = (_for_type*)ccv_get_dense_matrix_cell(cn, 0, 0, 0); \
43
0
  for (i = 0; i < rows * size; i++) \
44
0
  { \
45
0
    for (j = 0; j < cols * size; j++) \
46
0
    { \
47
0
      _for_type agv = agp[j * ch]; \
48
0
      _for_type mgv = mgp[j * ch]; \
49
0
      for (k = 1; k < ch; k++) \
50
0
        if (mgp[j * ch + k] > mgv) \
51
0
        { \
52
0
          mgv = mgp[j * ch + k]; \
53
0
          agv = agp[j * ch + k]; \
54
0
        } \
55
0
      _for_type agr0 = (ccv_clamp(agv, 0, 359.99) / 360.0) * (sbin * 2); \
56
0
      int ag0 = (int)agr0; \
57
0
      int ag1 = (ag0 + 1 < sbin * 2) ? ag0 + 1 : 0; \
58
0
      agr0 = agr0 - ag0; \
59
0
      _for_type agr1 = 1.0 - agr0; \
60
0
      mgv = mgv / 255.0; \
61
0
      _for_type yp = ((_for_type)i + 0.5) / (_for_type)size - 0.5; \
62
0
      _for_type xp = ((_for_type)j + 0.5) / (_for_type)size - 0.5; \
63
0
      int iyp = (int)floor(yp); \
64
0
      assert(iyp < rows); \
65
0
      int ixp = (int)floor(xp); \
66
0
      assert(ixp < cols); \
67
0
      _for_type vy0 = yp - iyp; \
68
0
      _for_type vx0 = xp - ixp; \
69
0
      _for_type vy1 = 1.0 - vy0; \
70
0
      _for_type vx1 = 1.0 - vx0; \
71
0
      if (ixp >= 0 && iyp >= 0) \
72
0
      { \
73
0
        cnp[iyp * cn->cols * sbin * 2 + ixp * sbin * 2 + ag0] += agr1 * vx1 * vy1 * mgv; \
74
0
        cnp[iyp * cn->cols * sbin * 2 + ixp * sbin * 2 + ag1] += agr0 * vx1 * vy1 * mgv; \
75
0
      } \
76
0
      if (ixp + 1 < cn->cols && iyp >= 0) \
77
0
      { \
78
0
        cnp[iyp * cn->cols * sbin * 2 + (ixp + 1) * sbin * 2 + ag0] += agr1 * vx0 * vy1 * mgv; \
79
0
        cnp[iyp * cn->cols * sbin * 2 + (ixp + 1) * sbin * 2 + ag1] += agr0 * vx0 * vy1 * mgv; \
80
0
      } \
81
0
      if (ixp >= 0 && iyp + 1 < cn->rows) \
82
0
      { \
83
0
        cnp[(iyp + 1) * cn->cols * sbin * 2 + ixp * sbin * 2 + ag0] += agr1 * vx1 * vy0 * mgv; \
84
0
        cnp[(iyp + 1) * cn->cols * sbin * 2 + ixp * sbin * 2 + ag1] += agr0 * vx1 * vy0 * mgv; \
85
0
      } \
86
0
      if (ixp + 1 < cn->cols && iyp + 1 < cn->rows) \
87
0
      { \
88
0
        cnp[(iyp + 1) * cn->cols * sbin * 2 + (ixp + 1) * sbin * 2 + ag0] += agr1 * vx0 * vy0 * mgv; \
89
0
        cnp[(iyp + 1) * cn->cols * sbin * 2 + (ixp + 1) * sbin * 2 + ag1] += agr0 * vx0 * vy0 * mgv; \
90
0
      } \
91
0
    } \
92
0
    agp += a->cols * ch; \
93
0
    mgp += a->cols * ch; \
94
0
  } \
95
0
  ccv_matrix_free(ag); \
96
0
  ccv_matrix_free(mg); \
97
0
  cnp = (_for_type*)ccv_get_dense_matrix_cell(cn, 0, 0, 0); \
98
0
  _for_type* cap = (_for_type*)ccv_get_dense_matrix_cell(ca, 0, 0, 0); \
99
0
  for (i = 0; i < rows; i++) \
100
0
  { \
101
0
    for (j = 0; j < cols; j++) \
102
0
    { \
103
0
      *cap = 0; \
104
0
      for (k = 0; k < sbin; k++) \
105
0
        *cap += (cnp[k] + cnp[k + sbin]) * (cnp[k] + cnp[k + sbin]); \
106
0
      cnp += 2 * sbin; \
107
0
      cap++; \
108
0
    } \
109
0
  } \
110
0
  cnp = (_for_type*)ccv_get_dense_matrix_cell(cn, 0, 0, 0); \
111
0
  cap = (_for_type*)ccv_get_dense_matrix_cell(ca, 0, 0, 0); \
112
0
  ccv_zero(db); \
113
0
  _for_type* dbp = (_for_type*)ccv_get_dense_matrix_cell(db, 0, 0, 0); \
114
0
  TNA(_for_type, 0, 1, cols + 1, cols, 0); \
115
0
  TNA(_for_type, 1, 1, 1, 0, 0); \
116
0
  TNA(_for_type, 2, 0, cols, cols, 0); \
117
0
  TNA(_for_type, 3, 0, 0, 0, 0); \
118
0
  cnp += 2 * sbin; \
119
0
  dbp += 3 * sbin + 4; \
120
0
  cap++; \
121
0
  for (j = 1; j < cols - 1; j++) \
122
0
  { \
123
0
    TNA(_for_type, 0, 1, cols + 1, cols, 0); \
124
0
    TNA(_for_type, 1, 1, 1, 0, 0); \
125
0
    TNA(_for_type, 2, -1, cols - 1, cols, 0); \
126
0
    TNA(_for_type, 3, -1, -1, 0, 0); \
127
0
    cnp += 2 * sbin; \
128
0
    dbp += 3 * sbin + 4; \
129
0
    cap++; \
130
0
  } \
131
0
  TNA(_for_type, 0, 0, cols, cols, 0); \
132
0
  TNA(_for_type, 1, 0, 0, 0, 0); \
133
0
  TNA(_for_type, 2, -1, cols - 1, cols, 0); \
134
0
  TNA(_for_type, 3, -1, -1, 0, 0); \
135
0
  cnp += 2 * sbin; \
136
0
  dbp += 3 * sbin + 4; \
137
0
  cap++; \
138
0
  for (i = 1; i < rows - 1; i++) \
139
0
  { \
140
0
    TNA(_for_type, 0, 1, cols + 1, cols, 0); \
141
0
    TNA(_for_type, 1, 1, -cols + 1, -cols, 0); \
142
0
    TNA(_for_type, 2, 0, cols, cols, 0); \
143
0
    TNA(_for_type, 3, 0, -cols, -cols, 0); \
144
0
    cnp += 2 * sbin; \
145
0
    dbp += 3 * sbin + 4; \
146
0
    cap++; \
147
0
    for (j = 1; j < cols - 1; j++) \
148
0
    { \
149
0
      TNA(_for_type, 0, 1, cols + 1, cols, 0); \
150
0
      TNA(_for_type, 1, 1, -cols + 1, -cols, 0); \
151
0
      TNA(_for_type, 2, -1, cols - 1, cols, 0); \
152
0
      TNA(_for_type, 3, -1, -cols - 1, -cols, 0); \
153
0
      cnp += 2 * sbin; \
154
0
      dbp += 3 * sbin + 4; \
155
0
      cap++; \
156
0
    } \
157
0
    TNA(_for_type, 0, 0, cols, cols, 0); \
158
0
    TNA(_for_type, 1, 0, -cols, -cols, 0); \
159
0
    TNA(_for_type, 2, -1, cols - 1, cols, 0); \
160
0
    TNA(_for_type, 3, -1, -cols - 1, -cols, 0); \
161
0
    cnp += 2 * sbin; \
162
0
    dbp += 3 * sbin + 4; \
163
0
    cap++; \
164
0
  } \
165
0
  TNA(_for_type, 0, 1, 1, 0, 0); \
166
0
  TNA(_for_type, 1, 1, -cols + 1, -cols, 0); \
167
0
  TNA(_for_type, 2, 0, 0, 0, 0); \
168
0
  TNA(_for_type, 3, 0, -cols, -cols, 0); \
169
0
  cnp += 2 * sbin; \
170
0
  dbp += 3 * sbin + 4; \
171
0
  cap++; \
172
0
  for (j = 1; j < cols - 1; j++) \
173
0
  { \
174
0
    TNA(_for_type, 0, 1, 1, 0, 0); \
175
0
    TNA(_for_type, 1, 1, -cols + 1, -cols, 0); \
176
0
    TNA(_for_type, 2, -1, -1, 0, 0); \
177
0
    TNA(_for_type, 3, -1, -cols - 1, -cols, 0); \
178
0
    cnp += 2 * sbin; \
179
0
    dbp += 3 * sbin + 4; \
180
0
    cap++; \
181
0
  } \
182
0
  TNA(_for_type, 0, 0, 0, 0, 0); \
183
0
  TNA(_for_type, 1, 0, -cols, -cols, 0); \
184
0
  TNA(_for_type, 2, -1, -1, 0, 0); \
185
0
  TNA(_for_type, 3, -1, -cols - 1, -cols, 0);
186
0
  ccv_matrix_typeof(db->type, for_block);
187
0
#undef for_block
188
0
#undef TNA
189
0
  ccv_matrix_free(cn);
190
0
  ccv_matrix_free(ca);
191
0
}
192
193
/* it is a supposely cleaner and faster implementation than original OpenCV (ccv_canny_deprecated,
194
 * removed, since the newer implementation achieve bit accuracy with OpenCV's), after a lot
195
 * profiling, the current implementation still uses integer to speed up */
196
void ccv_canny(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int size, double low_thresh, double high_thresh)
197
1
{
198
1
  assert(CCV_GET_CHANNEL(a->type) == CCV_C1);
199
1
  ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_canny(%d,%la,%la)", size, low_thresh, high_thresh), a->sig, CCV_EOF_SIGN);
200
1
  type = (type == 0) ? CCV_8U | CCV_C1 : 
CCV_GET_DATA_TYPE0
(type) | CCV_C10
;
201
1
  ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_C1 | CCV_ALL_DATA_TYPE, type, sig);
202
1
  ccv_object_return_if_cached(, db);
203
1
  if ((a->type & CCV_8U) || 
(a->type & CCV_32S)0
)
204
1
  {
205
1
    ccv_dense_matrix_t* dx = 0;
206
1
    ccv_dense_matrix_t* dy = 0;
207
1
    ccv_sobel(a, &dx, 0, size, 0);
208
1
    ccv_sobel(a, &dy, 0, 0, size);
209
1
    /* special case, all integer */
210
1
    int low = (int)(low_thresh + 0.5);
211
1
    int high = (int)(high_thresh + 0.5);
212
1
    int* dxi = dx->data.i32;
213
1
    int* dyi = dy->data.i32;
214
1
    int i, j;
215
1
    int* mbuf = (int*)alloca(3 * (a->cols + 2) * sizeof(int));
216
1
    memset(mbuf, 0, 3 * (a->cols + 2) * sizeof(int));
217
1
    int* rows[3];
218
1
    rows[0] = mbuf + 1;
219
1
    rows[1] = mbuf + (a->cols + 2) + 1;
220
1
    rows[2] = mbuf + 2 * (a->cols + 2) + 1;
221
501
    for (j = 0; j < a->cols; 
j++500
)
222
500
      rows[1][j] = abs(dxi[j]) + abs(dyi[j]);
223
1
    dxi += a->cols;
224
1
    dyi += a->cols;
225
1
    int* map = (int*)ccmalloc(sizeof(int) * (a->rows + 2) * (a->cols + 2));
226
1
    int map_cols = a->cols + 2;
227
1
    memset(map, 0, sizeof(int) * map_cols);
228
1
    int* map_ptr = map + map_cols + 1;
229
1
    int** stack = (int**)ccmalloc(sizeof(int*) * a->rows * a->cols);
230
1
    int** stack_top = stack;
231
1
    int** stack_bottom = stack;
232
501
    for (i = 1; i <= a->rows; 
i++500
)
233
500
    {
234
500
      /* the if clause should be unswitched automatically, no need to manually do so */
235
500
      if (i == a->rows)
236
1
        memset(rows[2], 0, sizeof(int) * a->cols);
237
499
      else
238
249k
        
for (j = 0; 499
j < a->cols;
j++249k
)
239
249k
          rows[2][j] = abs(dxi[j]) + abs(dyi[j]);
240
500
      int* _dx = dxi - a->cols;
241
500
      int* _dy = dyi - a->cols;
242
500
      map_ptr[-1] = 0;
243
500
      int suppress = 0;
244
250k
      for (j = 0; j < a->cols; 
j++250k
)
245
250k
      {
246
250k
        int f = rows[1][j];
247
250k
        if (f > low)
248
6.72k
        {
249
6.72k
          int x = abs(_dx[j]);
250
6.72k
          int y = abs(_dy[j]);
251
6.72k
          int s = _dx[j] ^ _dy[j];
252
6.72k
          /* x * tan(22.5) */
253
6.72k
          int tg22x = x * (int)(0.4142135623730950488016887242097 * (1 << 15) + 0.5);
254
6.72k
          /* x * tan(67.5) == 2 * x + x * tan(22.5) */
255
6.72k
          int tg67x = tg22x + ((x + x) << 15);
256
6.72k
          y <<= 15;
257
6.72k
          /* it is a little different from the Canny original paper because we adopted the coordinate system of
258
6.72k
           * top-left corner as origin. Thus, the derivative of y convolved with matrix:
259
6.72k
           * |-1 -2 -1|
260
6.72k
           * | 0  0  0|
261
6.72k
           * | 1  2  1|
262
6.72k
           * actually is the reverse of real y. Thus, the computed angle will be mirrored around x-axis.
263
6.72k
           * In this case, when angle is -45 (135), we compare with north-east and south-west, and for 45,
264
6.72k
           * we compare with north-west and south-east (in traditional coordinate system sense, the same if we
265
6.72k
           * adopt top-left corner as origin for "north", "south", "east", "west" accordingly) */
266
6.72k
#define high_block \
267
6.72k
          { \
268
1.14k
            if (f > high && !suppress && 
map_ptr[j - map_cols] != 2663
) \
269
1.14k
            { \
270
463
              map_ptr[j] = 2; \
271
463
              suppress = 1; \
272
463
              *(stack_top++) = map_ptr + j; \
273
682
            } else { \
274
682
              map_ptr[j] = 1; \
275
682
            } \
276
1.14k
            continue; \
277
1.14k
          }
278
6.72k
          /* sometimes, we end up with same f in integer domain, for that case, we will take the first occurrence
279
6.72k
           * suppressing the second with flag */
280
6.72k
          if (y < tg22x)
281
1.39k
          {
282
1.39k
            if (f > rows[1][j - 1] && 
f >= rows[1][j + 1]813
)
283
1.39k
              
high_block388
;
284
5.32k
          } else if (y > tg67x) {
285
1.34k
            if (f > rows[0][j] && 
f >= rows[2][j]786
)
286
1.34k
              
high_block357
;
287
3.98k
          } else {
288
3.98k
            s = s < 0 ? 
-12.45k
:
11.52k
;
289
3.98k
            if (f > rows[0][j - s] && 
f > rows[2][j + s]2.19k
)
290
3.98k
              
high_block400
;
291
3.58k
          }
292
6.72k
#undef high_block
293
6.72k
        }
294
250k
        map_ptr[j] = 0;
295
248k
        suppress = 0;
296
248k
      }
297
500
      map_ptr[a->cols] = 0;
298
500
      map_ptr += map_cols;
299
500
      dxi += a->cols;
300
500
      dyi += a->cols;
301
500
      int* row = rows[0];
302
500
      rows[0] = rows[1];
303
500
      rows[1] = rows[2];
304
500
      rows[2] = row;
305
500
    }
306
1
    memset(map_ptr - 1, 0, sizeof(int) * map_cols);
307
1
    int dr[] = {-1, 1, -map_cols - 1, -map_cols, -map_cols + 1, map_cols - 1, map_cols, map_cols + 1};
308
1.14k
    while (stack_top > stack_bottom)
309
1.14k
    {
310
1.14k
      map_ptr = *(--stack_top);
311
10.3k
      for (i = 0; i < 8; 
i++9.16k
)
312
9.16k
        if (map_ptr[dr[i]] == 1)
313
682
        {
314
682
          map_ptr[dr[i]] = 2;
315
682
          *(stack_top++) = map_ptr + dr[i];
316
682
        }
317
1.14k
    }
318
1
    map_ptr = map + map_cols + 1;
319
1
    unsigned char* b_ptr = db->data.u8;
320
1
#define for_block(_, _for_set) \
321
501
    
for (i = 0; 1
i < a->rows;
i++500
) \
322
500
    { \
323
250k
      for (j = 0; j < a->cols; 
j++250k
) \
324
500
        _for_set(b_ptr, j, (map_ptr[j] == 2)); \
325
500
      map_ptr += map_cols; \
326
500
      b_ptr += db->step; \
327
500
    }
328
1
    ccv_matrix_setter(db->type, for_block);
329
1
#undef for_block
330
1
    ccfree(stack);
331
1
    ccfree(map);
332
1
    ccv_matrix_free(dx);
333
1
    ccv_matrix_free(dy);
334
1
  } else {
335
0
    /* general case, use all ccv facilities to deal with it */
336
0
    ccv_dense_matrix_t* mg = 0;
337
0
    ccv_dense_matrix_t* ag = 0;
338
0
    ccv_gradient(a, &ag, 0, &mg, 0, size, size);
339
0
    ccv_matrix_free(ag);
340
0
    ccv_matrix_free(mg);
341
0
    /* FIXME: Canny implementation for general case */
342
0
  }
343
1
}
344
345
void ccv_close_outline(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type)
346
0
{
347
0
  assert((CCV_GET_CHANNEL(a->type) == CCV_C1) && ((a->type & CCV_8U) || (a->type & CCV_32S) || (a->type & CCV_64S)));
348
0
  ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_literal("ccv_close_outline"), a->sig, CCV_EOF_SIGN);
349
0
  type = ((type == 0) || (type & CCV_32F) || (type & CCV_64F)) ? CCV_GET_DATA_TYPE(a->type) | CCV_C1 : CCV_GET_DATA_TYPE(type) | CCV_C1;
350
0
  ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_C1 | CCV_ALL_DATA_TYPE, type, sig);
351
0
  ccv_object_return_if_cached(, db);
352
0
  int i, j;
353
0
  unsigned char* a_ptr = a->data.u8;
354
0
  unsigned char* b_ptr = db->data.u8;
355
0
  ccv_zero(db);
356
0
#define for_block(_for_get, _for_set_b, _for_get_b) \
357
0
  for (i = 0; i < a->rows - 1; i++) \
358
0
  { \
359
0
    for (j = 0; j < a->cols - 1; j++) \
360
0
    { \
361
0
      if (_for_get_b(b_ptr, j) == 0) \
362
0
        _for_set_b(b_ptr, j, _for_get(a_ptr, j)); \
363
0
      if (_for_get(a_ptr, j) != 0 && _for_get(a_ptr + a->step, j + 1) != 0) \
364
0
      { \
365
0
        _for_set_b(b_ptr + a->step, j, 1); \
366
0
        _for_set_b(b_ptr, j + 1, 1); \
367
0
      } \
368
0
      if (_for_get(a_ptr + a->step, j) != 0 && _for_get(a_ptr, j + 1) != 0) \
369
0
      { \
370
0
        _for_set_b(b_ptr, j, 1); \
371
0
        _for_set_b(b_ptr + a->step, j + 1, 1); \
372
0
      } \
373
0
    } \
374
0
    if (_for_get_b(b_ptr, a->cols - 1) == 0) \
375
0
      _for_set_b(b_ptr, a->cols - 1, _for_get(a_ptr, a->cols - 1)); \
376
0
    a_ptr += a->step; \
377
0
    b_ptr += db->step; \
378
0
  } \
379
0
  for (j = 0; j < a->cols; j++) \
380
0
  { \
381
0
    if (_for_get_b(b_ptr, j) == 0) \
382
0
      _for_set_b(b_ptr, j, _for_get(a_ptr, j)); \
383
0
  }
384
0
  ccv_matrix_getter_integer_only(a->type, ccv_matrix_setter_getter_integer_only, db->type, for_block);
385
0
#undef for_block
386
0
}
387
388
int ccv_otsu(ccv_dense_matrix_t* a, double* outvar, int range)
389
1
{
390
1
  assert((a->type & CCV_32S) || (a->type & CCV_8U));
391
1
  int* histogram = (int*)alloca(range * sizeof(int));
392
1
  memset(histogram, 0, sizeof(int) * range);
393
1
  int i, j;
394
1
  unsigned char* a_ptr = a->data.u8;
395
1
#define for_block(_, _for_get) \
396
7
  
for (i = 0; 1
i < a->rows;
i++6
) \
397
6
  { \
398
42
    for (j = 0; j < a->cols; 
j++36
) \
399
36
      histogram[ccv_clamp((int)_for_get(a_ptr, j), 0, range - 1)]++; \
400
6
    a_ptr += a->step; \
401
6
  }
402
1
  ccv_matrix_getter(a->type, for_block);
403
1
#undef for_block
404
1
  double sum = 0, sumB = 0;
405
7
  for (i = 0; i < range; 
i++6
)
406
6
    sum += i * histogram[i];
407
1
  int wB = 0, wF = 0, total = a->rows * a->cols;
408
1
  double maxVar = 0;
409
1
  int threshold = 0;
410
6
  for (i = 0; i < range; 
i++5
)
411
6
  {
412
6
    wB += histogram[i];
413
6
    if (wB == 0)
414
0
      continue;
415
6
    wF = total - wB;
416
6
    if (wF == 0)
417
1
      break;
418
5
    sumB += i * histogram[i];
419
5
    double mB = sumB / wB;
420
5
    double mF = (sum - sumB) / wF;
421
5
    double var = wB * wF * (mB - mF) * (mB - mF);
422
5
    if (var > maxVar)
423
3
    {
424
3
      maxVar = var;
425
3
      threshold = i;
426
3
    }
427
5
  }
428
1
  if (outvar != 0)
429
1
    *outvar = maxVar / total / total;
430
1
  return threshold;
431
1
}
432
433
0
#define LK_MAX_ITER (30)
434
0
#define LK_EPSILON (0.01)
435
436
/* this code is a rewrite from OpenCV's legendary Lucas-Kanade optical flow implementation */
437
void ccv_optical_flow_lucas_kanade(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_array_t* point_a, ccv_array_t** point_b, ccv_size_t win_size, int level, double min_eigen)
438
0
{
439
0
  assert(a && b && a->rows == b->rows && a->cols == b->cols);
440
0
  assert(CCV_GET_CHANNEL(a->type) == CCV_GET_CHANNEL(b->type) && CCV_GET_DATA_TYPE(a->type) == CCV_GET_DATA_TYPE(b->type));
441
0
  assert(CCV_GET_CHANNEL(a->type) == 1);
442
0
  assert(CCV_GET_DATA_TYPE(a->type) == CCV_8U);
443
0
  assert(point_a->rnum > 0);
444
0
  level = ccv_clamp(level + 1, 1, (int)(log((double)ccv_min(a->rows, a->cols) / ccv_max(win_size.width * 2, win_size.height * 2)) / log(2.0) + 0.5));
445
0
  ccv_declare_derived_signature(sig, a->sig != 0 && b->sig != 0 && point_a->sig != 0, ccv_sign_with_format(128, "ccv_optical_flow_lucas_kanade(%d,%d,%d,%la)", win_size.width, win_size.height, level, min_eigen), a->sig, b->sig, point_a->sig, CCV_EOF_SIGN);
446
0
  ccv_array_t* seq = *point_b = ccv_array_new(sizeof(ccv_decimal_point_with_status_t), point_a->rnum, sig);
447
0
  ccv_object_return_if_cached(, seq);
448
0
  seq->rnum = point_a->rnum;
449
0
  ccv_dense_matrix_t** pyr_a = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * level);
450
0
  ccv_dense_matrix_t** pyr_a_dx = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * level);
451
0
  ccv_dense_matrix_t** pyr_a_dy = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * level);
452
0
  ccv_dense_matrix_t** pyr_b = (ccv_dense_matrix_t**)alloca(sizeof(ccv_dense_matrix_t*) * level);
453
0
  int i, j, t, x, y;
454
0
  /* generating image pyramid */
455
0
  pyr_a[0] = a;
456
0
  pyr_a_dx[0] = pyr_a_dy[0] = 0;
457
0
  ccv_sobel(pyr_a[0], &pyr_a_dx[0], 0, 3, 0);
458
0
  ccv_sobel(pyr_a[0], &pyr_a_dy[0], 0, 0, 3);
459
0
  pyr_b[0] = b;
460
0
  for (i = 1; i < level; i++)
461
0
  {
462
0
    pyr_a[i] = pyr_a_dx[i] = pyr_a_dy[i] = pyr_b[i] = 0;
463
0
    ccv_sample_down(pyr_a[i - 1], &pyr_a[i], 0, 0, 0);
464
0
    ccv_sobel(pyr_a[i], &pyr_a_dx[i], 0, 3, 0);
465
0
    ccv_sobel(pyr_a[i], &pyr_a_dy[i], 0, 0, 3);
466
0
    ccv_sample_down(pyr_b[i - 1], &pyr_b[i], 0, 0, 0);
467
0
  }
468
0
  int* wi = (int*)alloca(sizeof(int) * win_size.width * win_size.height);
469
0
  int* widx = (int*)alloca(sizeof(int) * win_size.width * win_size.height);
470
0
  int* widy = (int*)alloca(sizeof(int) * win_size.width * win_size.height);
471
0
  ccv_decimal_point_t half_win = ccv_decimal_point((win_size.width - 1) * 0.5f, (win_size.height - 1) * 0.5f);
472
0
  const int W_BITS14 = 14, W_BITS7 = 7, W_BITS9 = 9;
473
0
  const float FLT_SCALE = 1.0f / (1 << 25);
474
0
  // clean up status to 1
475
0
  for (i = 0; i < point_a->rnum; i++)
476
0
  {
477
0
    ccv_decimal_point_with_status_t* point_with_status = (ccv_decimal_point_with_status_t*)ccv_array_get(seq, i);
478
0
    point_with_status->status = 1;
479
0
  }
480
0
  int prev_rows, prev_cols;
481
0
  for (t = level - 1; t >= 0; t--)
482
0
  {
483
0
    ccv_dense_matrix_t* a = pyr_a[t];
484
0
    ccv_dense_matrix_t* adx = pyr_a_dx[t];
485
0
    ccv_dense_matrix_t* ady = pyr_a_dy[t];
486
0
    assert(CCV_GET_DATA_TYPE(adx->type) == CCV_32S);
487
0
    assert(CCV_GET_DATA_TYPE(ady->type) == CCV_32S);
488
0
    ccv_dense_matrix_t* b = pyr_b[t];
489
0
    for (i = 0; i < point_a->rnum; i++)
490
0
    {
491
0
      ccv_decimal_point_t prev_point = *(ccv_decimal_point_t*)ccv_array_get(point_a, i);
492
0
      ccv_decimal_point_with_status_t* point_with_status = (ccv_decimal_point_with_status_t*)ccv_array_get(seq, i);
493
0
      prev_point.x = prev_point.x / (float)(1 << t);
494
0
      prev_point.y = prev_point.y / (float)(1 << t);
495
0
      ccv_decimal_point_t next_point;
496
0
      if (t == level - 1)
497
0
        next_point = prev_point;
498
0
      else {
499
0
        next_point.x = point_with_status->point.x * 2 + (a->cols - prev_cols * 2) * 0.5;
500
0
        next_point.y = point_with_status->point.y * 2 + (a->rows - prev_rows * 2) * 0.5;
501
0
      }
502
0
      point_with_status->point = next_point;
503
0
      prev_point.x -= half_win.x;
504
0
      prev_point.y -= half_win.y;
505
0
      ccv_point_t iprev_point = ccv_point((int)prev_point.x, (int)prev_point.y);
506
0
      if (iprev_point.x < 0 || iprev_point.x >= a->cols - win_size.width - 1 ||
507
0
        iprev_point.y < 0 || iprev_point.y >= a->rows - win_size.height - 1)
508
0
      {
509
0
        if (t == 0)
510
0
          point_with_status->status = 0;
511
0
        continue;
512
0
      }
513
0
      float xd = prev_point.x - iprev_point.x;
514
0
      float yd = prev_point.y - iprev_point.y;
515
0
      int iw00 = (int)((1 - xd) * (1 - yd) * (1 << W_BITS14) + 0.5);
516
0
      int iw01 = (int)(xd * (1 - yd) * (1 << W_BITS14) + 0.5);
517
0
      int iw10 = (int)((1 - xd) * yd * (1 << W_BITS14) + 0.5);
518
0
      int iw11 = (1 << W_BITS14) - iw00 - iw01 - iw10;
519
0
      float a11 = 0, a12 = 0, a22 = 0;
520
0
      unsigned char* a_ptr = (unsigned char*)ccv_get_dense_matrix_cell_by(CCV_C1 | CCV_8U, a, iprev_point.y, iprev_point.x, 0);
521
0
      int* adx_ptr = (int*)ccv_get_dense_matrix_cell_by(CCV_C1 | CCV_32S, adx, iprev_point.y, iprev_point.x, 0);
522
0
      int* ady_ptr = (int*)ccv_get_dense_matrix_cell_by(CCV_C1 | CCV_32S, ady, iprev_point.y, iprev_point.x, 0);
523
0
      int* wi_ptr = wi;
524
0
      int* widx_ptr = widx;
525
0
      int* widy_ptr = widy;
526
0
      for (y = 0; y < win_size.height; y++)
527
0
      {
528
0
        for (x = 0; x < win_size.width; x++)
529
0
        {
530
0
          wi_ptr[x] = ccv_descale(a_ptr[x] * iw00 + a_ptr[x + 1] * iw01 + a_ptr[x + a->step] * iw10 + a_ptr[x + a->step + 1] * iw11, W_BITS7);
531
0
          // because we use 3x3 sobel, which scaled derivative up by 4
532
0
          widx_ptr[x] = ccv_descale(adx_ptr[x] * iw00 + adx_ptr[x + 1] * iw01 + adx_ptr[x + adx->cols] * iw10 + adx_ptr[x + adx->cols + 1] * iw11, W_BITS9);
533
0
          widy_ptr[x] = ccv_descale(ady_ptr[x] * iw00 + ady_ptr[x + 1] * iw01 + ady_ptr[x + ady->cols] + iw10 + ady_ptr[x + ady->cols + 1] * iw11, W_BITS9);
534
0
          a11 += (float)(widx_ptr[x] * widx_ptr[x]);
535
0
          a12 += (float)(widx_ptr[x] * widy_ptr[x]);
536
0
          a22 += (float)(widy_ptr[x] * widy_ptr[x]);
537
0
        }
538
0
        a_ptr += a->step;
539
0
        adx_ptr += adx->cols;
540
0
        ady_ptr += ady->cols;
541
0
        wi_ptr += win_size.width;
542
0
        widx_ptr += win_size.width;
543
0
        widy_ptr += win_size.width;
544
0
      }
545
0
      a11 *= FLT_SCALE;
546
0
      a12 *= FLT_SCALE;
547
0
      a22 *= FLT_SCALE;
548
0
      float D = a11 * a22 - a12 * a12;
549
0
      float eigen = (a22 + a11 - sqrtf((a11 - a22) * (a11 - a22) + 4.0f * a12 * a12)) / (2 * win_size.width * win_size.height);
550
0
      if (eigen < min_eigen || D < FLT_EPSILON)
551
0
      {
552
0
        if (t == 0)
553
0
          point_with_status->status = 0;
554
0
        continue;
555
0
      }
556
0
      D = 1.0f / D;
557
0
      next_point.x -= half_win.x;
558
0
      next_point.y -= half_win.y;
559
0
      ccv_decimal_point_t prev_delta;
560
0
      for (j = 0; j < LK_MAX_ITER; j++)
561
0
      {
562
0
        ccv_point_t inext_point = ccv_point((int)next_point.x, (int)next_point.y);
563
0
        if (inext_point.x < 0 || inext_point.x >= a->cols - win_size.width - 1 ||
564
0
          inext_point.y < 0 || inext_point.y >= a->rows - win_size.height - 1)
565
0
          break;
566
0
        float xd = next_point.x - inext_point.x;
567
0
        float yd = next_point.y - inext_point.y;
568
0
        int iw00 = (int)((1 - xd) * (1 - yd) * (1 << W_BITS14) + 0.5);
569
0
        int iw01 = (int)(xd * (1 - yd) * (1 << W_BITS14) + 0.5);
570
0
        int iw10 = (int)((1 - xd) * yd * (1 << W_BITS14) + 0.5);
571
0
        int iw11 = (1 << W_BITS14) - iw00 - iw01 - iw10;
572
0
        float b1 = 0, b2 = 0;
573
0
        unsigned char* b_ptr = (unsigned char*)ccv_get_dense_matrix_cell_by(CCV_C1 | CCV_8U, b, inext_point.y, inext_point.x, 0);
574
0
        int* wi_ptr = wi;
575
0
        int* widx_ptr = widx;
576
0
        int* widy_ptr = widy;
577
0
        for (y = 0; y < win_size.height; y++)
578
0
        {
579
0
          for (x = 0; x < win_size.width; x++)
580
0
          {
581
0
            int diff = ccv_descale(b_ptr[x] * iw00 + b_ptr[x + 1] * iw01 + b_ptr[x + b->step] * iw10 + b_ptr[x + b->step + 1] * iw11, W_BITS7) - wi_ptr[x];
582
0
            b1 += (float)(diff * widx_ptr[x]);
583
0
            b2 += (float)(diff * widy_ptr[x]);
584
0
          }
585
0
          b_ptr += b->step;
586
0
          wi_ptr += win_size.width;
587
0
          widx_ptr += win_size.width;
588
0
          widy_ptr += win_size.width;
589
0
        }
590
0
        b1 *= FLT_SCALE;
591
0
        b2 *= FLT_SCALE;
592
0
        ccv_decimal_point_t delta = ccv_decimal_point((a12 * b2 - a22 * b1) * D, (a12 * b1 - a11 * b2) * D);
593
0
        next_point.x += delta.x;
594
0
        next_point.y += delta.y;
595
0
        if (delta.x * delta.x + delta.y * delta.y < LK_EPSILON)
596
0
          break;
597
0
        if (j > 0 && fabs(prev_delta.x - delta.x) < 0.01 && fabs(prev_delta.y - delta.y) < 0.01)
598
0
        {
599
0
          next_point.x -= delta.x * 0.5;
600
0
          next_point.y -= delta.y * 0.5;
601
0
          break;
602
0
        }
603
0
        prev_delta = delta;
604
0
      }
605
0
      ccv_point_t inext_point = ccv_point((int)next_point.x, (int)next_point.y);
606
0
      if (inext_point.x < 0 || inext_point.x >= a->cols - win_size.width - 1 ||
607
0
        inext_point.y < 0 || inext_point.y >= a->rows - win_size.height - 1)
608
0
        point_with_status->status = 0;
609
0
      else {
610
0
        point_with_status->point.x = next_point.x + half_win.x;
611
0
        point_with_status->point.y = next_point.y + half_win.y;
612
0
      }
613
0
    }
614
0
    prev_rows = a->rows;
615
0
    prev_cols = a->cols;
616
0
    ccv_matrix_free(adx);
617
0
    ccv_matrix_free(ady);
618
0
    if (t > 0)
619
0
    {
620
0
      ccv_matrix_free(a);
621
0
      ccv_matrix_free(b);
622
0
    }
623
0
  }
624
0
}