Bug Summary

File:ccv_numeric.c
Warning:line 1283, column 33
The left operand of '+' is a garbage value

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name ccv_numeric.c -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -pic-is-pie -mframe-pointer=none -menable-no-infs -menable-no-nans -fapprox-func -funsafe-math-optimizations -fno-signed-zeros -mreassociate -freciprocal-math -ffp-contract=fast -fno-rounding-math -ffast-math -ffinite-math-only -complex-range=basic -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -target-feature +sse2 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/home/liu/actions-runner/_work/ccv/ccv/lib -fcoverage-compilation-dir=/home/liu/actions-runner/_work/ccv/ccv/lib -resource-dir /usr/local/lib/clang/19 -I . -I /usr/local/cuda/include -D HAVE_CBLAS -D HAVE_LIBPNG -D HAVE_LIBJPEG -D HAVE_FFTW3 -D HAVE_PTHREAD -D HAVE_LIBLINEAR -D HAVE_TESSERACT -D HAVE_AVCODEC -D HAVE_AVFORMAT -D HAVE_AVUTIL -D HAVE_SWSCALE -D HAVE_SSE2 -D HAVE_GSL -D HAVE_CUDA -D HAVE_CUDNN -D HAVE_NCCL -D USE_SYSTEM_CUB -I /usr/local/include -internal-isystem /usr/local/lib/clang/19/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/12/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -O3 -ferror-limit 19 -fgnuc-version=4.2.1 -fskip-odr-check-in-gmf -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /home/liu/actions-runner/_work/ccv/ccv/_analyze/2024-11-14-215813-628724-1 -x c ccv_numeric.c
1#include "ccv.h"
2#include "ccv_internal.h"
3#include <complex_Complex.h>
4#ifdef HAVE_FFTW31
5#include <pthread.h>
6#include <fftw3.h>
7#else
8#include "3rdparty/kissfft/kiss_fftndr.h"
9#include "3rdparty/kissfft/kissf_fftndr.h"
10#endif
11
12const ccv_minimize_param_t ccv_minimize_default_params = {
13 .interp = 0.1,
14 .extrap = 3.0,
15 .max_iter = 20,
16 .ratio = 10.0,
17 .sig = 0.1,
18 .rho = 0.05,
19};
20
21void ccv_invert(ccv_matrix_t* a, ccv_matrix_t** b, int type)
22{
23}
24
25void ccv_solve(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type)
26{
27}
28
29void ccv_eigen(ccv_dense_matrix_t* a, ccv_dense_matrix_t** vector, ccv_dense_matrix_t** lambda, int type, double epsilon)
30{
31 ccv_declare_derived_signature(vsig, a->sig != 0, ccv_sign_with_literal("ccv_eigen(vector)"), a->sig, CCV_EOF_SIGN)char _ccv_identifier_31[] = ("ccv_eigen(vector)"); size_t _ccv_string_size_31
= sizeof(_ccv_identifier_31);; uint64_t vsig = (a->sig !=
0) ? ccv_cache_generate_signature(_ccv_identifier_31, _ccv_string_size_31
, a->sig, ((uint64_t)0)) : 0;;
32 ccv_declare_derived_signature(lsig, a->sig != 0, ccv_sign_with_literal("ccv_eigen(lambda)"), a->sig, CCV_EOF_SIGN)char _ccv_identifier_32[] = ("ccv_eigen(lambda)"); size_t _ccv_string_size_32
= sizeof(_ccv_identifier_32);; uint64_t lsig = (a->sig !=
0) ? ccv_cache_generate_signature(_ccv_identifier_32, _ccv_string_size_32
, a->sig, ((uint64_t)0)) : 0;;
33 assert(CCV_GET_CHANNEL(a->type) == 1)((void) sizeof ((((a->type) & 0xFFF) == 1) ? 1 : 0), __extension__
({ if (((a->type) & 0xFFF) == 1) ; else __assert_fail
("CCV_GET_CHANNEL(a->type) == 1", "ccv_numeric.c", 33, __extension__
__PRETTY_FUNCTION__); }));
34 type = (type == 0) ? CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) | CCV_C1 : CCV_GET_DATA_TYPE(type)((type) & 0xFF000) | CCV_C1;
35 // as of now, this function only support real symmetric matrix
36 ccv_dense_matrix_t* dvector = *vector = ccv_dense_matrix_renew(*vector, a->rows, a->cols, CCV_32F | CCV_64F | CCV_C1, type, vsig);
37 ccv_dense_matrix_t* dlambda = *lambda = ccv_dense_matrix_renew(*lambda, 1, a->cols, CCV_32F | CCV_64F | CCV_C1, type, lsig);
38 assert(CCV_GET_DATA_TYPE(dvector->type) == CCV_GET_DATA_TYPE(dlambda->type))((void) sizeof ((((dvector->type) & 0xFF000) == ((dlambda
->type) & 0xFF000)) ? 1 : 0), __extension__ ({ if (((dvector
->type) & 0xFF000) == ((dlambda->type) & 0xFF000
)) ; else __assert_fail ("CCV_GET_DATA_TYPE(dvector->type) == CCV_GET_DATA_TYPE(dlambda->type)"
, "ccv_numeric.c", 38, __extension__ __PRETTY_FUNCTION__); })
);
39 ccv_object_return_if_cached(, dvector, dlambda){ if ((!(dvector) || (((int*)(dvector))[0] & CCV_GARBAGE)
) && (!(dlambda) || (((int*)(dlambda))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE))) { (void)((dvector) && (
((int*)(dvector))[0] &= ~CCV_GARBAGE));(void)((dlambda) &&
(((int*)(dlambda))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));;; return ; } };
40 double* ja = (double*)ccmallocmalloc(sizeof(double) * a->rows * a->cols);
41 int i, j;
42 unsigned char* aptr = a->data.u8;
43 assert(a->rows > 0 && a->cols > 0)((void) sizeof ((a->rows > 0 && a->cols >
0) ? 1 : 0), __extension__ ({ if (a->rows > 0 &&
a->cols > 0) ; else __assert_fail ("a->rows > 0 && a->cols > 0"
, "ccv_numeric.c", 43, __extension__ __PRETTY_FUNCTION__); })
);
44#define for_block(_, _for_get) \
45 for (i = 0; i < a->rows; i++) \
46 { \
47 for (j = 0; j < a->cols; j++) \
48 ja[i * a->cols + j] = _for_get(aptr, j); \
49 aptr += a->step; \
50 }
51 ccv_matrix_getter(a->type, for_block){ switch (((a->type) & 0xFF000)) { case CCV_32S: { for_block
(, _ccv_get_32s_value); break; } case CCV_32F: { for_block(, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(, _ccv_get_64f_value); break; }
default: { for_block(, _ccv_get_8u_value); } } };
52#undef for_block
53 ccv_zero(dvector);
54 ccv_zero(dlambda);
55 unsigned char* dvptr = dvector->data.u8;
56#define for_block(_, _for_set) \
57 for (i = 0; i < a->cols; i++) \
58 _for_set(dvptr, i * a->cols + i, 1);
59 ccv_matrix_setter(dvector->type, for_block){ switch (((dvector->type) & 0xFF000)) { case CCV_32S:
{ for_block(, _ccv_set_32s_value); break; } case CCV_32F: { for_block
(, _ccv_set_32f_value); break; } case CCV_64S: { for_block(, _ccv_set_64s_value
); break; } case CCV_64F: { for_block(, _ccv_set_64f_value); break
; } default: { for_block(, _ccv_set_8u_value); } } };
60#undef for_block
61 double accuracy = 0;
62 for (i = 0; i < a->rows * a->cols; i++)
63 accuracy += ja[i];
64 accuracy = sqrt(2 * accuracy);
65 int p, q;
66 unsigned char* dlptr = dlambda->data.u8;
67 int flag = 1;
68 assert(a->rows == a->cols)((void) sizeof ((a->rows == a->cols) ? 1 : 0), __extension__
({ if (a->rows == a->cols) ; else __assert_fail ("a->rows == a->cols"
, "ccv_numeric.c", 68, __extension__ __PRETTY_FUNCTION__); })
);
69#define for_block(_, _for_set, _for_get) \
70 do { \
71 if (!flag) \
72 accuracy = accuracy * 0.5; \
73 flag = 0; \
74 for (p = 0; p < a->rows; p++) \
75 { \
76 for (q = p + 1; q < a->cols; q++) \
77 if (fabs(ja[p * a->cols + q]) > accuracy) \
78 { \
79 double x = -ja[p * a->cols + q]; \
80 double y = (ja[q * a->cols + q] - ja[p * a->cols + p]) * 0.5; \
81 double omega = (x == 0 && y == 0) ? 1 : x / sqrt(x * x + y * y); \
82 if (y < 0) \
83 omega = -omega; \
84 double sn = 1.0 + sqrt(1.0 - omega * omega); \
85 sn = omega / sqrt(2 * sn); \
86 double cn = sqrt(1.0 - sn * sn); \
87 double fpp = ja[p * a->cols + p]; \
88 double fpq = ja[p * a->cols + q]; \
89 double fqq = ja[q * a->cols + q]; \
90 ja[p * a->cols + p] = fpp * cn * cn + fqq * sn * sn + fpq * omega; \
91 ja[q * a->cols + q] = fpp * sn * sn + fqq * cn * cn - fpq * omega; \
92 ja[p * a->cols + q] = ja[q * a->cols + p] = 0; \
93 for (i = 0; i < a->cols; i++) \
94 if (i != q && i != p) \
95 { \
96 fpp = ja[p * a->cols + i]; \
97 fqq = ja[q * a->cols + i]; \
98 ja[p * a->cols + i] = fpp * cn + fqq * sn; \
99 ja[q * a->cols + i] = -fpp * sn + fqq * cn; \
100 } \
101 for (i = 0; i < a->rows; i++) \
102 if (i != q && i != p) \
103 { \
104 fpp = ja[i * a->cols + p]; \
105 fqq = ja[i * a->cols + q]; \
106 ja[i * a->cols + p] = fpp * cn + fqq * sn; \
107 ja[i * a->cols + q] = -fpp * sn + fqq * cn; \
108 } \
109 for (i = 0; i < a->cols; i++) \
110 { \
111 fpp = _for_get(dvptr, p * a->cols + i); \
112 fqq = _for_get(dvptr, q * a->cols + i); \
113 _for_set(dvptr, p * a->cols + i, fpp * cn + fqq * sn); \
114 _for_set(dvptr, q * a->cols + i, -fpp * sn + fqq * cn); \
115 } \
116 for (i = 0; i < a->cols; i++) \
117 _for_set(dlptr, i, ja[i * a->cols + i]); \
118 flag = 1; \
119 break; \
120 } \
121 if (flag) \
122 break; \
123 } \
124 } while (accuracy > epsilon);
125 ccv_matrix_setter_getter(dvector->type, for_block){ switch (((dvector->type) & 0xFF000)) { case CCV_32S:
{ for_block(, _ccv_set_32s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(, _ccv_set_32f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(, _ccv_set_64s_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(, _ccv_set_64f_value, _ccv_get_64f_value
); break; } default: { for_block(, _ccv_set_8u_value, _ccv_get_8u_value
); } } };
126#undef for_block
127 ccfreefree(ja);
128}
129
130void ccv_minimize(ccv_dense_matrix_t* x, int length, double red, ccv_minimize_f func, ccv_minimize_param_t params, void* data)
131{
132 ccv_dense_matrix_t* df0 = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
133 ccv_zero(df0);
134 ccv_dense_matrix_t* df3 = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
135 ccv_zero(df3);
136 ccv_dense_matrix_t* dF0 = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
137 ccv_zero(dF0);
138 ccv_dense_matrix_t* s = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
139 ccv_zero(s);
140 ccv_dense_matrix_t* x0 = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
141 ccv_zero(x0);
142 ccv_dense_matrix_t* xn = ccv_dense_matrix_new(x->rows, x->cols, x->type, 0, 0);
143 ccv_zero(xn);
144
145 double F0 = 0, f0 = 0, f1 = 0, f2 = 0, f3 = 0, f4 = 0;
146 double x1 = 0, x2 = 0, x3 = 0, x4 = 0;
147 double d0 = 0, d1 = 0, d2 = 0, d3 = 0, d4 = 0;
148 double A = 0, B = 0;
149 int ls_failed = 0;
150
151 int i, j, k;
152 func(x, &f0, df0, data);
153 d0 = 0;
154 unsigned char* df0p = df0->data.u8;
155 unsigned char* sp = s->data.u8;
156 for (i = 0; i < x->rows; i++)
157 {
158 for (j = 0; j < x->cols; j++)
159 {
160 double ss = ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)]))));
161 ccv_set_value(x->type, sp, j, -ss, 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(sp))[(j)] = (int)(-ss) >> 0; break; case CCV_32F: ((
float*)(sp))[(j)] = (float)-ss; break; case CCV_64S: ((int64_t
*)(sp))[(j)] = (int64_t)(-ss) >> 0; break; case CCV_64F
: ((double*)(sp))[(j)] = (double)-ss; break; default: ((unsigned
char*)(sp))[(j)] = ({ typeof (0) _a = (0); typeof (255) _b =
(255); typeof ((int)(-ss) >> 0) _x = ((int)(-ss) >>
0); (_x < _a) ? _a : ((_x > _b) ? _b : _x); }); };
162 d0 += -ss * ss;
163 }
164 df0p += x->step;
165 sp += x->step;
166 }
167 x3 = red / (1.0 - d0);
168 int l = (length > 0) ? length : -length;
169 int ls = (length > 0) ? 1 : 0;
170 int eh = (length > 0) ? 0 : 1;
171 for (k = 0; k < l;)
172 {
173 k += ls;
174 memcpy(x0->data.u8, x->data.u8, x->rows * x->step);
175 memcpy(dF0->data.u8, df0->data.u8, x->rows * x->step);
176 F0 = f0;
177 int m = ccv_min(params.max_iter, (length > 0) ? params.max_iter : l - k)({ typeof (params.max_iter) _a = (params.max_iter); typeof ((
length > 0) ? params.max_iter : l - k) _b = ((length > 0
) ? params.max_iter : l - k); (_a < _b) ? _a : _b; });
178 for (;;)
179 {
180 x2 = 0;
181 f3 = f2 = f0;
182 d2 = d0;
183 memcpy(df3->data.u8, df0->data.u8, x->rows * x->step);
184 while (m > 0)
185 {
186 m--;
187 k += eh;
188 unsigned char* sp = s->data.u8;
189 unsigned char* xp = x->data.u8;
190 unsigned char* xnp = xn->data.u8;
191 for (i = 0; i < x->rows; i++)
192 {
193 for (j = 0; j < x->cols; j++)
194 ccv_set_value(x->type, xnp, j, x3 * ccv_get_value(x->type, sp, j) + ccv_get_value(x->type, xp, j), 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(xnp))[(j)] = (int)(x3 * (((x->type) & CCV_32S) ? ((
int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*)(sp
))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp))[(j
)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)] : (
(unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0; break; case CCV_32F
: ((float*)(xnp))[(j)] = (float)x3 * (((x->type) & CCV_32S
) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)
] : ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)])))); break; case CCV_64S: ((int64_t
*)(xnp))[(j)] = (int64_t)(x3 * (((x->type) & CCV_32S) ?
((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*
)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp)
)[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)]
: ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0; break; case CCV_64F
: ((double*)(xnp))[(j)] = (double)x3 * (((x->type) & CCV_32S
) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)
] : ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)])))); break; default: ((unsigned
char*)(xnp))[(j)] = ({ typeof (0) _a = (0); typeof (255) _b =
(255); typeof ((int)(x3 * (((x->type) & CCV_32S) ? ((
int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*)(sp
))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp))[(j
)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)] : (
(unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0) _x = ((int)(
x3 * (((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x
->type) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type
) & CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)
])))) + (((x->type) & CCV_32S) ? ((int*)(xp))[(j)] : (
((x->type) & CCV_32F) ? ((float*)(xp))[(j)] : (((x->
type) & CCV_64S) ? ((int64_t*)(xp))[(j)] : (((x->type)
& CCV_64F) ? ((double*)(xp))[(j)] : ((unsigned char*)(xp
))[(j)]))))) >> 0); (_x < _a) ? _a : ((_x > _b) ?
_b : _x); }); };
195 sp += x->step;
196 xp += x->step;
197 xnp += x->step;
198 }
199 if (func(xn, &f3, df3, data) == 0)
200 break;
201 else
202 x3 = (x2 + x3) * 0.5;
203 }
204 if (f3 < F0)
205 {
206 memcpy(x0->data.u8, xn->data.u8, x->rows * x->step);
207 memcpy(dF0->data.u8, df3->data.u8, x->rows * x->step);
208 F0 = f3;
209 }
210 d3 = 0;
211 unsigned char* df3p = df3->data.u8;
212 unsigned char* sp = s->data.u8;
213 for (i = 0; i < x->rows; i++)
214 {
215 for (j = 0; j < x->cols; j++)
216 d3 += ccv_get_value(x->type, df3p, j)(((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)])))) * ccv_get_value(x->type, sp, j)(((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type) &
CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) & CCV_64F
) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)]))));
217 df3p += x->step;
218 sp += x->step;
219 }
220 if ((d3 > params.sig * d0) || (f3 > f0 + x3 * params.rho * d0) || (m <= 0))
221 break;
222 x1 = x2; f1 = f2; d1 = d2;
223 x2 = x3; f2 = f3; d2 = d3;
224 A = 6.0 * (f1 - f2) + 3.0 * (d2 + d1) * (x2 - x1);
225 B = 3.0 * (f2 - f1) - (2.0 * d1 + d2) * (x2 - x1);
226 x3 = B * B - A * d1 * (x2 - x1);
227 if (x3 < 0)
228 x3 = x2 * params.extrap;
229 else {
230 x3 = x1 - d1 * (x2 - x1) * (x2 - x1) / (B + sqrt(x3));
231 if (x3 < 0)
232 x3 = x2 * params.extrap;
233 else {
234 if (x3 > x2 * params.extrap)
235 x3 = x2 * params.extrap;
236 else if (x3 < x2 + params.interp * (x2 - x1))
237 x3 = x2 + params.interp * (x2 - x1);
238 }
239 }
240 }
241 while (((fabs(d3) > -params.sig * d0) || (f3 > f0 + x3 * params.rho * d0)) && (m > 0))
242 {
243 if ((d3 > 1e-8) || (f3 > f0 + x3 * params.rho * d0))
244 {
245 x4 = x3; f4 = f3; d4 = d3;
246 } else {
247 x2 = x3; f2 = f3; d2 = d3;
248 }
249 if (f4 > f0)
250 x3 = x2 - (0.5 * d2 * (x4 - x2) * (x4 - x2)) / (f4 - f2 - d2 * (x4 - x2));
251 else {
252 A = 6.0 * (f2 - f4) / (x4 - x2) + 3.0 * (d4 + d2);
253 B = 3.0 * (f4 - f2) - (2.0 * d2 + d4) * (x4 - x2);
254 x3 = B * B - A * d2 * (x4 - x2) * (x4 - x2);
255 x3 = (x3 < 0) ? (x2 + x4) * 0.5 : x2 + (sqrt(x3) - B) / A;
256 }
257 x3 = ccv_max(ccv_min(x3, x4 - params.interp * (x4 - x2)), x2 + params.interp * (x4 - x2))({ typeof (({ typeof (x3) _a = (x3); typeof (x4 - params.interp
* (x4 - x2)) _b = (x4 - params.interp * (x4 - x2)); (_a <
_b) ? _a : _b; })) _a = (({ typeof (x3) _a = (x3); typeof (x4
- params.interp * (x4 - x2)) _b = (x4 - params.interp * (x4 -
x2)); (_a < _b) ? _a : _b; })); typeof (x2 + params.interp
* (x4 - x2)) _b = (x2 + params.interp * (x4 - x2)); (_a >
_b) ? _a : _b; });
258 sp = s->data.u8;
259 unsigned char* xp = x->data.u8;
260 unsigned char* xnp = xn->data.u8;
261 for (i = 0; i < x->rows; i++)
262 {
263 for (j = 0; j < x->cols; j++)
264 ccv_set_value(x->type, xnp, j, x3 * ccv_get_value(x->type, sp, j) + ccv_get_value(x->type, xp, j), 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(xnp))[(j)] = (int)(x3 * (((x->type) & CCV_32S) ? ((
int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*)(sp
))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp))[(j
)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)] : (
(unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0; break; case CCV_32F
: ((float*)(xnp))[(j)] = (float)x3 * (((x->type) & CCV_32S
) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)
] : ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)])))); break; case CCV_64S: ((int64_t
*)(xnp))[(j)] = (int64_t)(x3 * (((x->type) & CCV_32S) ?
((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*
)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp)
)[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)]
: ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0; break; case CCV_64F
: ((double*)(xnp))[(j)] = (double)x3 * (((x->type) & CCV_32S
) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)
] : ((unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)])))); break; default: ((unsigned
char*)(xnp))[(j)] = ({ typeof (0) _a = (0); typeof (255) _b =
(255); typeof ((int)(x3 * (((x->type) & CCV_32S) ? ((
int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*)(sp
))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp))[(j
)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)] : (
(unsigned char*)(sp))[(j)])))) + (((x->type) & CCV_32S
) ? ((int*)(xp))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(xp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(xp
))[(j)] : (((x->type) & CCV_64F) ? ((double*)(xp))[(j)
] : ((unsigned char*)(xp))[(j)]))))) >> 0) _x = ((int)(
x3 * (((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x
->type) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type
) & CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)
])))) + (((x->type) & CCV_32S) ? ((int*)(xp))[(j)] : (
((x->type) & CCV_32F) ? ((float*)(xp))[(j)] : (((x->
type) & CCV_64S) ? ((int64_t*)(xp))[(j)] : (((x->type)
& CCV_64F) ? ((double*)(xp))[(j)] : ((unsigned char*)(xp
))[(j)]))))) >> 0); (_x < _a) ? _a : ((_x > _b) ?
_b : _x); }); };
265 sp += x->step;
266 xp += x->step;
267 xnp += x->step;
268 }
269 func(xn, &f3, df3, data);
270 if (f3 < F0)
271 {
272 memcpy(x0->data.u8, xn->data.u8, x->rows * x->step);
273 memcpy(dF0->data.u8, df3->data.u8, x->rows * x->step);
274 F0 = f3;
275 }
276 m--;
277 k += eh;
278 d3 = 0;
279 sp = s->data.u8;
280 unsigned char* df3p = df3->data.u8;
281 for (i = 0; i < x->rows; i++)
282 {
283 for (j = 0; j < x->cols; j++)
284 d3 += ccv_get_value(x->type, df3p, j)(((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)])))) * ccv_get_value(x->type, sp, j)(((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type) &
CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) & CCV_64F
) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)]))));
285 df3p += x->step;
286 sp += x->step;
287 }
288 }
289 if ((fabs(d3) < -params.sig * d0) && (f3 < f0 + x3 * params.rho * d0))
290 {
291 memcpy(x->data.u8, xn->data.u8, x->rows * x->step);
292 f0 = f3;
293 double df0_df3 = 0;
294 double df3_df3 = 0;
295 double df0_df0 = 0;
296 unsigned char* df0p = df0->data.u8;
297 unsigned char* df3p = df3->data.u8;
298 for (i = 0; i < x->rows; i++)
299 {
300 for (j = 0; j < x->cols; j++)
301 {
302 df0_df0 += ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)])))) * ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)]))));
303 df0_df3 += ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)])))) * ccv_get_value(x->type, df3p, j)(((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)]))));
304 df3_df3 += ccv_get_value(x->type, df3p, j)(((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)])))) * ccv_get_value(x->type, df3p, j)(((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)]))));
305 }
306 df0p += x->step;
307 df3p += x->step;
308 }
309 double slr = (df3_df3 - df0_df3) / df0_df0;
310 df3p = df3->data.u8;
311 unsigned char* sp = s->data.u8;
312 for (i = 0; i < x->rows; i++)
313 {
314 for (j = 0; j < x->cols; j++)
315 ccv_set_value(x->type, sp, j, slr * ccv_get_value(x->type, sp, j) - ccv_get_value(x->type, df3p, j), 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(sp))[(j)] = (int)(slr * (((x->type) & CCV_32S) ? ((
int*)(sp))[(j)] : (((x->type) & CCV_32F) ? ((float*)(sp
))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(sp))[(j
)] : (((x->type) & CCV_64F) ? ((double*)(sp))[(j)] : (
(unsigned char*)(sp))[(j)])))) - (((x->type) & CCV_32S
) ? ((int*)(df3p))[(j)] : (((x->type) & CCV_32F) ? ((float
*)(df3p))[(j)] : (((x->type) & CCV_64S) ? ((int64_t*)(
df3p))[(j)] : (((x->type) & CCV_64F) ? ((double*)(df3p
))[(j)] : ((unsigned char*)(df3p))[(j)]))))) >> 0; break
; case CCV_32F: ((float*)(sp))[(j)] = (float)slr * (((x->type
) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F
) ? ((float*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t
*)(sp))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp)
)[(j)] : ((unsigned char*)(sp))[(j)])))) - (((x->type) &
CCV_32S) ? ((int*)(df3p))[(j)] : (((x->type) & CCV_32F
) ? ((float*)(df3p))[(j)] : (((x->type) & CCV_64S) ? (
(int64_t*)(df3p))[(j)] : (((x->type) & CCV_64F) ? ((double
*)(df3p))[(j)] : ((unsigned char*)(df3p))[(j)])))); break; case
CCV_64S: ((int64_t*)(sp))[(j)] = (int64_t)(slr * (((x->type
) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F
) ? ((float*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t
*)(sp))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp)
)[(j)] : ((unsigned char*)(sp))[(j)])))) - (((x->type) &
CCV_32S) ? ((int*)(df3p))[(j)] : (((x->type) & CCV_32F
) ? ((float*)(df3p))[(j)] : (((x->type) & CCV_64S) ? (
(int64_t*)(df3p))[(j)] : (((x->type) & CCV_64F) ? ((double
*)(df3p))[(j)] : ((unsigned char*)(df3p))[(j)]))))) >> 0
; break; case CCV_64F: ((double*)(sp))[(j)] = (double)slr * (
((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->type
) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type) &
CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) & CCV_64F
) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)])))) - (
((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)])))); break; default: ((unsigned char*)(sp))[(j)] = ({ typeof
(0) _a = (0); typeof (255) _b = (255); typeof ((int)(slr * (
((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->type
) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type) &
CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) & CCV_64F
) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)])))) - (
((x->type) & CCV_32S) ? ((int*)(df3p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df3p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df3p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df3p))[(j)] : ((unsigned char*)(df3p))
[(j)]))))) >> 0) _x = ((int)(slr * (((x->type) &
CCV_32S) ? ((int*)(sp))[(j)] : (((x->type) & CCV_32F)
? ((float*)(sp))[(j)] : (((x->type) & CCV_64S) ? ((int64_t
*)(sp))[(j)] : (((x->type) & CCV_64F) ? ((double*)(sp)
)[(j)] : ((unsigned char*)(sp))[(j)])))) - (((x->type) &
CCV_32S) ? ((int*)(df3p))[(j)] : (((x->type) & CCV_32F
) ? ((float*)(df3p))[(j)] : (((x->type) & CCV_64S) ? (
(int64_t*)(df3p))[(j)] : (((x->type) & CCV_64F) ? ((double
*)(df3p))[(j)] : ((unsigned char*)(df3p))[(j)]))))) >> 0
); (_x < _a) ? _a : ((_x > _b) ? _b : _x); }); };
316 df3p += x->step;
317 sp += x->step;
318 }
319 memcpy(df0->data.u8, df3->data.u8, x->rows * x->step);
320 d3 = d0;
321 d0 = 0;
322 df0p = df0->data.u8;
323 sp = s->data.u8;
324 for (i = 0; i < x->rows; i++)
325 {
326 for (j = 0; j < x->cols; j++)
327 {
328 d0 += ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)])))) * ccv_get_value(x->type, sp, j)(((x->type) & CCV_32S) ? ((int*)(sp))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(sp))[(j)] : (((x->type) &
CCV_64S) ? ((int64_t*)(sp))[(j)] : (((x->type) & CCV_64F
) ? ((double*)(sp))[(j)] : ((unsigned char*)(sp))[(j)]))));
329 }
330 df0p += x->step;
331 sp += x->step;
332 }
333 if (d0 > 0)
334 {
335 d0 = 0;
336 df0p = df0->data.u8;
337 sp = s->data.u8;
338 for (i = 0; i < x->rows; i++)
339 {
340 for (j = 0; j < x->cols; j++)
341 {
342 double ss = ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)]))));
343 ccv_set_value(x->type, sp, j, -ss, 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(sp))[(j)] = (int)(-ss) >> 0; break; case CCV_32F: ((
float*)(sp))[(j)] = (float)-ss; break; case CCV_64S: ((int64_t
*)(sp))[(j)] = (int64_t)(-ss) >> 0; break; case CCV_64F
: ((double*)(sp))[(j)] = (double)-ss; break; default: ((unsigned
char*)(sp))[(j)] = ({ typeof (0) _a = (0); typeof (255) _b =
(255); typeof ((int)(-ss) >> 0) _x = ((int)(-ss) >>
0); (_x < _a) ? _a : ((_x > _b) ? _b : _x); }); };
344 d0 += -ss * ss;
345 }
346 df0p += x->step;
347 sp += x->step;
348 }
349 }
350 x3 = x3 * ccv_min(params.ratio, d3 / (d0 - 1e-8))({ typeof (params.ratio) _a = (params.ratio); typeof (d3 / (d0
- 1e-8)) _b = (d3 / (d0 - 1e-8)); (_a < _b) ? _a : _b; });
351 ls_failed = 0;
352 } else {
353 memcpy(x->data.u8, x0->data.u8, x->rows * x->step);
354 memcpy(df0->data.u8, dF0->data.u8, x->rows * x->step);
355 f0 = F0;
356 if (ls_failed)
357 break;
358 d0 = 0;
359 unsigned char* df0p = df0->data.u8;
360 unsigned char* sp = s->data.u8;
361 for (i = 0; i < x->rows; i++)
362 {
363 for (j = 0; j < x->cols; j++)
364 {
365 double ss = ccv_get_value(x->type, df0p, j)(((x->type) & CCV_32S) ? ((int*)(df0p))[(j)] : (((x->
type) & CCV_32F) ? ((float*)(df0p))[(j)] : (((x->type)
& CCV_64S) ? ((int64_t*)(df0p))[(j)] : (((x->type) &
CCV_64F) ? ((double*)(df0p))[(j)] : ((unsigned char*)(df0p))
[(j)]))));
366 ccv_set_value(x->type, sp, j, -ss, 0)switch ((((x->type)) & 0xFF000)) { case CCV_32S: ((int
*)(sp))[(j)] = (int)(-ss) >> 0; break; case CCV_32F: ((
float*)(sp))[(j)] = (float)-ss; break; case CCV_64S: ((int64_t
*)(sp))[(j)] = (int64_t)(-ss) >> 0; break; case CCV_64F
: ((double*)(sp))[(j)] = (double)-ss; break; default: ((unsigned
char*)(sp))[(j)] = ({ typeof (0) _a = (0); typeof (255) _b =
(255); typeof ((int)(-ss) >> 0) _x = ((int)(-ss) >>
0); (_x < _a) ? _a : ((_x > _b) ? _b : _x); }); };
367 d0 += -ss * ss;
368 }
369 df0p += x->step;
370 sp += x->step;
371 }
372 x3 = red / (1.0 - d0);
373 ls_failed = 1;
374 }
375 }
376 ccv_matrix_free(s);
377 ccv_matrix_free(x0);
378 ccv_matrix_free(xn);
379 ccv_matrix_free(dF0);
380 ccv_matrix_free(df0);
381 ccv_matrix_free(df3);
382}
383
384#ifdef HAVE_FFTW31
385/* optimal FFT size table is adopted from OpenCV */
386static const int _ccv_optimal_fft_size[] = {
387 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
388 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160,
389 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375,
390 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720,
391 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,
392 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875,
393 1920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592,
394 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
395 3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400,
396 5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290,
397 7500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000,
398 10125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500,
399 12800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000,
400 16200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683,
401 20000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576,
402 25000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720,
403 31104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864,
404 37500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080,
405 46656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675,
406 55296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800,
407 65536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125,
408 78732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312,
409 93750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000,
410 109350, 110592, 112500, 115200, 116640, 118098, 120000, 121500, 122880, 124416, 125000,
411 128000, 129600, 131072, 131220, 135000, 138240, 139968, 140625, 144000, 145800, 147456,
412 150000, 151875, 153600, 155520, 156250, 157464, 160000, 162000, 163840, 164025, 165888,
413 168750, 172800, 174960, 177147, 180000, 182250, 184320, 186624, 187500, 192000, 194400,
414 196608, 196830, 200000, 202500, 204800, 207360, 209952, 216000, 218700, 221184, 225000,
415 230400, 233280, 234375, 236196, 240000, 243000, 245760, 248832, 250000, 253125, 256000,
416 259200, 262144, 262440, 270000, 273375, 276480, 279936, 281250, 288000, 291600, 294912,
417 295245, 300000, 303750, 307200, 311040, 312500, 314928, 320000, 324000, 327680, 328050,
418 331776, 337500, 345600, 349920, 354294, 360000, 364500, 368640, 373248, 375000, 384000,
419 388800, 390625, 393216, 393660, 400000, 405000, 409600, 414720, 419904, 421875, 432000,
420 437400, 442368, 450000, 455625, 460800, 466560, 468750, 472392, 480000, 486000, 491520,
421 492075, 497664, 500000, 506250, 512000, 518400, 524288, 524880, 531441, 540000, 546750,
422 552960, 559872, 562500, 576000, 583200, 589824, 590490, 600000, 607500, 614400, 622080,
423 625000, 629856, 640000, 648000, 655360, 656100, 663552, 675000, 691200, 699840, 703125,
424 708588, 720000, 729000, 737280, 746496, 750000, 759375, 768000, 777600, 781250, 786432,
425 787320, 800000, 810000, 819200, 820125, 829440, 839808, 843750, 864000, 874800, 884736,
426 885735, 900000, 911250, 921600, 933120, 937500, 944784, 960000, 972000, 983040, 984150,
427 995328, 1000000, 1012500, 1024000, 1036800, 1048576, 1049760, 1062882, 1080000, 1093500,
428 1105920, 1119744, 1125000, 1152000, 1166400, 1171875, 1179648, 1180980, 1200000,
429 1215000, 1228800, 1244160, 1250000, 1259712, 1265625, 1280000, 1296000, 1310720,
430 1312200, 1327104, 1350000, 1366875, 1382400, 1399680, 1406250, 1417176, 1440000,
431 1458000, 1474560, 1476225, 1492992, 1500000, 1518750, 1536000, 1555200, 1562500,
432 1572864, 1574640, 1594323, 1600000, 1620000, 1638400, 1640250, 1658880, 1679616,
433 1687500, 1728000, 1749600, 1769472, 1771470, 1800000, 1822500, 1843200, 1866240,
434 1875000, 1889568, 1920000, 1944000, 1953125, 1966080, 1968300, 1990656, 2000000,
435 2025000, 2048000, 2073600, 2097152, 2099520, 2109375, 2125764, 2160000, 2187000,
436 2211840, 2239488, 2250000, 2278125, 2304000, 2332800, 2343750, 2359296, 2361960,
437 2400000, 2430000, 2457600, 2460375, 2488320, 2500000, 2519424, 2531250, 2560000,
438 2592000, 2621440, 2624400, 2654208, 2657205, 2700000, 2733750, 2764800, 2799360,
439 2812500, 2834352, 2880000, 2916000, 2949120, 2952450, 2985984, 3000000, 3037500,
440 3072000, 3110400, 3125000, 3145728, 3149280, 3188646, 3200000, 3240000, 3276800,
441 3280500, 3317760, 3359232, 3375000, 3456000, 3499200, 3515625, 3538944, 3542940,
442 3600000, 3645000, 3686400, 3732480, 3750000, 3779136, 3796875, 3840000, 3888000,
443 3906250, 3932160, 3936600, 3981312, 4000000, 4050000, 4096000, 4100625, 4147200,
444 4194304, 4199040, 4218750, 4251528, 4320000, 4374000, 4423680, 4428675, 4478976,
445 4500000, 4556250, 4608000, 4665600, 4687500, 4718592, 4723920, 4782969, 4800000,
446 4860000, 4915200, 4920750, 4976640, 5000000, 5038848, 5062500, 5120000, 5184000,
447 5242880, 5248800, 5308416, 5314410, 5400000, 5467500, 5529600, 5598720, 5625000,
448 5668704, 5760000, 5832000, 5859375, 5898240, 5904900, 5971968, 6000000, 6075000,
449 6144000, 6220800, 6250000, 6291456, 6298560, 6328125, 6377292, 6400000, 6480000,
450 6553600, 6561000, 6635520, 6718464, 6750000, 6834375, 6912000, 6998400, 7031250,
451 7077888, 7085880, 7200000, 7290000, 7372800, 7381125, 7464960, 7500000, 7558272,
452 7593750, 7680000, 7776000, 7812500, 7864320, 7873200, 7962624, 7971615, 8000000,
453 8100000, 8192000, 8201250, 8294400, 8388608, 8398080, 8437500, 8503056, 8640000,
454 8748000, 8847360, 8857350, 8957952, 9000000, 9112500, 9216000, 9331200, 9375000,
455 9437184, 9447840, 9565938, 9600000, 9720000, 9765625, 9830400, 9841500, 9953280,
456 10000000, 10077696, 10125000, 10240000, 10368000, 10485760, 10497600, 10546875, 10616832,
457 10628820, 10800000, 10935000, 11059200, 11197440, 11250000, 11337408, 11390625, 11520000,
458 11664000, 11718750, 11796480, 11809800, 11943936, 12000000, 12150000, 12288000, 12301875,
459 12441600, 12500000, 12582912, 12597120, 12656250, 12754584, 12800000, 12960000, 13107200,
460 13122000, 13271040, 13286025, 13436928, 13500000, 13668750, 13824000, 13996800, 14062500,
461 14155776, 14171760, 14400000, 14580000, 14745600, 14762250, 14929920, 15000000, 15116544,
462 15187500, 15360000, 15552000, 15625000, 15728640, 15746400, 15925248, 15943230, 16000000,
463 16200000, 16384000, 16402500, 16588800, 16777216, 16796160, 16875000, 17006112, 17280000,
464 17496000, 17578125, 17694720, 17714700, 17915904, 18000000, 18225000, 18432000, 18662400,
465 18750000, 18874368, 18895680, 18984375, 19131876, 19200000, 19440000, 19531250, 19660800,
466 19683000, 19906560, 20000000, 20155392, 20250000, 20480000, 20503125, 20736000, 20971520,
467 20995200, 21093750, 21233664, 21257640, 21600000, 21870000, 22118400, 22143375, 22394880,
468 22500000, 22674816, 22781250, 23040000, 23328000, 23437500, 23592960, 23619600, 23887872,
469 23914845, 24000000, 24300000, 24576000, 24603750, 24883200, 25000000, 25165824, 25194240,
470 25312500, 25509168, 25600000, 25920000, 26214400, 26244000, 26542080, 26572050, 26873856,
471 27000000, 27337500, 27648000, 27993600, 28125000, 28311552, 28343520, 28800000, 29160000,
472 29296875, 29491200, 29524500, 29859840, 30000000, 30233088, 30375000, 30720000, 31104000,
473 31250000, 31457280, 31492800, 31640625, 31850496, 31886460, 32000000, 32400000, 32768000,
474 32805000, 33177600, 33554432, 33592320, 33750000, 34012224, 34171875, 34560000, 34992000,
475 35156250, 35389440, 35429400, 35831808, 36000000, 36450000, 36864000, 36905625, 37324800,
476 37500000, 37748736, 37791360, 37968750, 38263752, 38400000, 38880000, 39062500, 39321600,
477 39366000, 39813120, 39858075, 40000000, 40310784, 40500000, 40960000, 41006250, 41472000,
478 41943040, 41990400, 42187500, 42467328, 42515280, 43200000, 43740000, 44236800, 44286750,
479 44789760, 45000000, 45349632, 45562500, 46080000, 46656000, 46875000, 47185920, 47239200,
480 47775744, 47829690, 48000000, 48600000, 48828125, 49152000, 49207500, 49766400, 50000000,
481 50331648, 50388480, 50625000, 51018336, 51200000, 51840000, 52428800, 52488000, 52734375,
482 53084160, 53144100, 53747712, 54000000, 54675000, 55296000, 55987200, 56250000, 56623104,
483 56687040, 56953125, 57600000, 58320000, 58593750, 58982400, 59049000, 59719680, 60000000,
484 60466176, 60750000, 61440000, 61509375, 62208000, 62500000, 62914560, 62985600, 63281250,
485 63700992, 63772920, 64000000, 64800000, 65536000, 65610000, 66355200, 66430125, 67108864,
486 67184640, 67500000, 68024448, 68343750, 69120000, 69984000, 70312500, 70778880, 70858800,
487 71663616, 72000000, 72900000, 73728000, 73811250, 74649600, 75000000, 75497472, 75582720,
488 75937500, 76527504, 76800000, 77760000, 78125000, 78643200, 78732000, 79626240, 79716150,
489 80000000, 80621568, 81000000, 81920000, 82012500, 82944000, 83886080, 83980800, 84375000,
490 84934656, 85030560, 86400000, 87480000, 87890625, 88473600, 88573500, 89579520, 90000000,
491 90699264, 91125000, 92160000, 93312000, 93750000, 94371840, 94478400, 94921875, 95551488,
492 95659380, 96000000, 97200000, 97656250, 98304000, 98415000, 99532800, 100000000,
493 100663296, 100776960, 101250000, 102036672, 102400000, 102515625, 103680000, 104857600,
494 104976000, 105468750, 106168320, 106288200, 107495424, 108000000, 109350000, 110592000,
495 110716875, 111974400, 112500000, 113246208, 113374080, 113906250, 115200000, 116640000,
496 117187500, 117964800, 118098000, 119439360, 119574225, 120000000, 120932352, 121500000,
497 122880000, 123018750, 124416000, 125000000, 125829120, 125971200, 126562500, 127401984,
498 127545840, 128000000, 129600000, 131072000, 131220000, 132710400, 132860250, 134217728,
499 134369280, 135000000, 136048896, 136687500, 138240000, 139968000, 140625000, 141557760,
500 141717600, 143327232, 144000000, 145800000, 146484375, 147456000, 147622500, 149299200,
501 150000000, 150994944, 151165440, 151875000, 153055008, 153600000, 155520000, 156250000,
502 157286400, 157464000, 158203125, 159252480, 159432300, 160000000, 161243136, 162000000,
503 163840000, 164025000, 165888000, 167772160, 167961600, 168750000, 169869312, 170061120,
504 170859375, 172800000, 174960000, 175781250, 176947200, 177147000, 179159040, 180000000,
505 181398528, 182250000, 184320000, 184528125, 186624000, 187500000, 188743680, 188956800,
506 189843750, 191102976, 191318760, 192000000, 194400000, 195312500, 196608000, 196830000,
507 199065600, 199290375, 200000000, 201326592, 201553920, 202500000, 204073344, 204800000,
508 205031250, 207360000, 209715200, 209952000, 210937500, 212336640, 212576400, 214990848,
509 216000000, 218700000, 221184000, 221433750, 223948800, 225000000, 226492416, 226748160,
510 227812500, 230400000, 233280000, 234375000, 235929600, 236196000, 238878720, 239148450,
511 240000000, 241864704, 243000000, 244140625, 245760000, 246037500, 248832000, 250000000,
512 251658240, 251942400, 253125000, 254803968, 255091680, 256000000, 259200000, 262144000,
513 262440000, 263671875, 265420800, 265720500, 268435456, 268738560, 270000000, 272097792,
514 273375000, 276480000, 279936000, 281250000, 283115520, 283435200, 284765625, 286654464,
515 288000000, 291600000, 292968750, 294912000, 295245000, 298598400, 300000000, 301989888,
516 302330880, 303750000, 306110016, 307200000, 307546875, 311040000, 312500000, 314572800,
517 314928000, 316406250, 318504960, 318864600, 320000000, 322486272, 324000000, 327680000,
518 328050000, 331776000, 332150625, 335544320, 335923200, 337500000, 339738624, 340122240,
519 341718750, 345600000, 349920000, 351562500, 353894400, 354294000, 358318080, 360000000,
520 362797056, 364500000, 368640000, 369056250, 373248000, 375000000, 377487360, 377913600,
521 379687500, 382205952, 382637520, 384000000, 388800000, 390625000, 393216000, 393660000,
522 398131200, 398580750, 400000000, 402653184, 403107840, 405000000, 408146688, 409600000,
523 410062500, 414720000, 419430400, 419904000, 421875000, 424673280, 425152800, 429981696,
524 432000000, 437400000, 439453125, 442368000, 442867500, 447897600, 450000000, 452984832,
525 453496320, 455625000, 460800000, 466560000, 468750000, 471859200, 472392000, 474609375,
526 477757440, 478296900, 480000000, 483729408, 486000000, 488281250, 491520000, 492075000,
527 497664000, 500000000, 503316480, 503884800, 506250000, 509607936, 510183360, 512000000,
528 512578125, 518400000, 524288000, 524880000, 527343750, 530841600, 531441000, 536870912,
529 537477120, 540000000, 544195584, 546750000, 552960000, 553584375, 559872000, 562500000,
530 566231040, 566870400, 569531250, 573308928, 576000000, 583200000, 585937500, 589824000,
531 590490000, 597196800, 597871125, 600000000, 603979776, 604661760, 607500000, 612220032,
532 614400000, 615093750, 622080000, 625000000, 629145600, 629856000, 632812500, 637009920,
533 637729200, 640000000, 644972544, 648000000, 655360000, 656100000, 663552000, 664301250,
534 671088640, 671846400, 675000000, 679477248, 680244480, 683437500, 691200000, 699840000,
535 703125000, 707788800, 708588000, 716636160, 720000000, 725594112, 729000000, 732421875,
536 737280000, 738112500, 746496000, 750000000, 754974720, 755827200, 759375000, 764411904,
537 765275040, 768000000, 777600000, 781250000, 786432000, 787320000, 791015625, 796262400,
538 797161500, 800000000, 805306368, 806215680, 810000000, 816293376, 819200000, 820125000,
539 829440000, 838860800, 839808000, 843750000, 849346560, 850305600, 854296875, 859963392,
540 864000000, 874800000, 878906250, 884736000, 885735000, 895795200, 900000000, 905969664,
541 906992640, 911250000, 921600000, 922640625, 933120000, 937500000, 943718400, 944784000,
542 949218750, 955514880, 956593800, 960000000, 967458816, 972000000, 976562500, 983040000,
543 984150000, 995328000, 996451875, 1000000000, 1006632960, 1007769600, 1012500000,
544 1019215872, 1020366720, 1024000000, 1025156250, 1036800000, 1048576000, 1049760000,
545 1054687500, 1061683200, 1062882000, 1073741824, 1074954240, 1080000000, 1088391168,
546 1093500000, 1105920000, 1107168750, 1119744000, 1125000000, 1132462080, 1133740800,
547 1139062500, 1146617856, 1152000000, 1166400000, 1171875000, 1179648000, 1180980000,
548 1194393600, 1195742250, 1200000000, 1207959552, 1209323520, 1215000000, 1220703125,
549 1224440064, 1228800000, 1230187500, 1244160000, 1250000000, 1258291200, 1259712000,
550 1265625000, 1274019840, 1275458400, 1280000000, 1289945088, 1296000000, 1310720000,
551 1312200000, 1318359375, 1327104000, 1328602500, 1342177280, 1343692800, 1350000000,
552 1358954496, 1360488960, 1366875000, 1382400000, 1399680000, 1406250000, 1415577600,
553 1417176000, 1423828125, 1433272320, 1440000000, 1451188224, 1458000000, 1464843750,
554 1474560000, 1476225000, 1492992000, 1500000000, 1509949440, 1511654400, 1518750000,
555 1528823808, 1530550080, 1536000000, 1537734375, 1555200000, 1562500000, 1572864000,
556 1574640000, 1582031250, 1592524800, 1594323000, 1600000000, 1610612736, 1612431360,
557 1620000000, 1632586752, 1638400000, 1640250000, 1658880000, 1660753125, 1677721600,
558 1679616000, 1687500000, 1698693120, 1700611200, 1708593750, 1719926784, 1728000000,
559 1749600000, 1757812500, 1769472000, 1771470000, 1791590400, 1800000000, 1811939328,
560 1813985280, 1822500000, 1843200000, 1845281250, 1866240000, 1875000000, 1887436800,
561 1889568000, 1898437500, 1911029760, 1913187600, 1920000000, 1934917632, 1944000000,
562 1953125000, 1966080000, 1968300000, 1990656000, 1992903750, 2000000000, 2013265920,
563 2015539200, 2025000000, 2038431744, 2040733440, 2048000000, 2050312500, 2073600000,
564 2097152000, 2099520000, 2109375000, 2123366400, 2125764000
565};
566
567static int _ccv_get_optimal_fft_size(int size)
568{
569 int a = 0, b = sizeof(_ccv_optimal_fft_size) / sizeof(_ccv_optimal_fft_size[0]) - 1;
570 if((unsigned)size >= (unsigned)_ccv_optimal_fft_size[b])
571 return -1;
572 while(a < b)
573 {
574 int c = (a + b) >> 1;
575 if(size <= _ccv_optimal_fft_size[c])
576 b = c;
577 else
578 a = c + 1;
579 }
580 return _ccv_optimal_fft_size[b];
581}
582
583static pthread_mutex_t fftw_plan_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, PTHREAD_MUTEX_TIMED_NP, 0, 0, { 0, 0 } } };
584
585static void _ccv_filter_fftw(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t* d, int padding_pattern)
586{
587 int ch = CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF);
588 int fft_type = (CCV_GET_DATA_TYPE(d->type)((d->type) & 0xFF000) == CCV_8U || CCV_GET_DATA_TYPE(d->type)((d->type) & 0xFF000) == CCV_32F) ? CCV_32F : CCV_64F;
589 int rows = ccv_min(a->rows + b->rows - 1, _ccv_get_optimal_fft_size(b->rows * 3))({ typeof (a->rows + b->rows - 1) _a = (a->rows + b->
rows - 1); typeof (_ccv_get_optimal_fft_size(b->rows * 3))
_b = (_ccv_get_optimal_fft_size(b->rows * 3)); (_a < _b
) ? _a : _b; });
590 int cols = ccv_min(a->cols + b->cols - 1, _ccv_get_optimal_fft_size(b->cols * 3))({ typeof (a->cols + b->cols - 1) _a = (a->cols + b->
cols - 1); typeof (_ccv_get_optimal_fft_size(b->cols * 3))
_b = (_ccv_get_optimal_fft_size(b->cols * 3)); (_a < _b
) ? _a : _b; });
591 int cols_2c = 2 * (cols / 2 + 1);
592 void* fftw_a;
593 void* fftw_b;
594 void* fftw_d;
595 fftw_plan p, pinv;
596 fftwf_plan pf, pinvf;
597 if (fft_type == CCV_32F)
598 {
599 fftw_a = fftwf_malloc(rows * cols_2c * ch * sizeof(float));
600 fftw_b = fftwf_malloc(rows * cols_2c * ch * sizeof(float));
601 fftw_d = fftwf_malloc(rows * cols_2c * ch * sizeof(float));
602 pthread_mutex_lock(&fftw_plan_mutex);
603 if (ch == 1)
604 {
605 pf = fftwf_plan_dft_r2c_2d(rows, cols, 0, 0, FFTW_ESTIMATE(1U << 6));
606 pinvf = fftwf_plan_dft_c2r_2d(rows, cols, 0, 0, FFTW_ESTIMATE(1U << 6));
607 } else {
608 int ndim[] = {rows, cols};
609 pf = fftwf_plan_many_dft_r2c(2, ndim, ch, 0, 0, ch, 1, 0, 0, ch, 1, FFTW_ESTIMATE(1U << 6));
610 pinvf = fftwf_plan_many_dft_c2r(2, ndim, ch, 0, 0, ch, 1, 0, 0, ch, 1, FFTW_ESTIMATE(1U << 6));
611 }
612 pthread_mutex_unlock(&fftw_plan_mutex);
613 } else {
614 fftw_a = fftw_malloc(rows * cols_2c * ch * sizeof(double));
615 fftw_b = fftw_malloc(rows * cols_2c * ch * sizeof(double));
616 fftw_d = fftw_malloc(rows * cols_2c * ch * sizeof(double));
617 pthread_mutex_lock(&fftw_plan_mutex);
618 if (ch == 1)
619 {
620 p = fftw_plan_dft_r2c_2d(rows, cols, 0, 0, FFTW_ESTIMATE(1U << 6));
621 pinv = fftw_plan_dft_c2r_2d(rows, cols, 0, 0, FFTW_ESTIMATE(1U << 6));
622 } else {
623 int ndim[] = {rows, cols};
624 p = fftw_plan_many_dft_r2c(2, ndim, ch, 0, 0, ch, 1, 0, 0, ch, 1, FFTW_ESTIMATE(1U << 6));
625 pinv = fftw_plan_many_dft_c2r(2, ndim, ch, 0, 0, ch, 1, 0, 0, ch, 1, FFTW_ESTIMATE(1U << 6));
626 }
627 pthread_mutex_unlock(&fftw_plan_mutex);
628 }
629 memset(fftw_b, 0, rows * cols_2c * ch * CCV_GET_DATA_TYPE_SIZE(fft_type)_ccv_get_data_type_size[((fft_type) & 0xFF000) >> 12
]);
630
631 int i, j, k;
632 unsigned char* m_ptr = b->data.u8;
633 // to flip matrix b is crucial, this problem only shows when I changed to a more sophisticated test case
634#define for_block(_for_type, _for_get) \
635 _for_type* fftw_ptr = (_for_type*)fftw_b + (b->rows - 1) * cols_2c * ch; \
636 for (i = 0; i < b->rows; i++) \
637 { \
638 for (j = 0; j < b->cols; j++) \
639 for (k = 0; k < ch; k++) \
640 fftw_ptr[(b->cols - 1 - j) * ch + k] = _for_get(m_ptr, j * ch + k); \
641 fftw_ptr -= cols_2c * ch; \
642 m_ptr += b->step; \
643 }
644 ccv_matrix_typeof(fft_type, ccv_matrix_getter, b->type, for_block){ switch (((fft_type) & 0xFF000)) { case CCV_32S: { { switch
(((b->type) & 0xFF000)) { case CCV_32S: { for_block(int
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(int
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(int
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(int
, _ccv_get_64f_value); break; } default: { for_block(int, _ccv_get_8u_value
); } } }; break; } case CCV_32F: { { switch (((b->type) &
0xFF000)) { case CCV_32S: { for_block(float, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(float, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(float, _ccv_get_64f_value
); break; } default: { for_block(float, _ccv_get_8u_value); }
} }; break; } case CCV_64S: { { switch (((b->type) & 0xFF000
)) { case CCV_32S: { for_block(int64_t, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(int64_t, _ccv_get_32f_value); break
; } case CCV_64S: { for_block(int64_t, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(int64_t, _ccv_get_64f_value); break
; } default: { for_block(int64_t, _ccv_get_8u_value); } } }; break
; } case CCV_64F: { { switch (((b->type) & 0xFF000)) {
case CCV_32S: { for_block(double, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(double, _ccv_get_32f_value); break
; } case CCV_64S: { for_block(double, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, _ccv_get_64f_value); break
; } default: { for_block(double, _ccv_get_8u_value); } } }; break
; } default: { { switch (((b->type) & 0xFF000)) { case
CCV_32S: { for_block(unsigned char, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(unsigned char, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(unsigned char, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(unsigned char, _ccv_get_64f_value
); break; } default: { for_block(unsigned char, _ccv_get_8u_value
); } } }; } } };
645#undef for_block
646 if (fft_type == CCV_32F)
647 fftwf_execute_dft_r2c(pf, (float*)fftw_b, (fftwf_complex*)fftw_b);
648 else
649 fftw_execute_dft_r2c(p, (double*)fftw_b, (fftw_complex*)fftw_b);
650 /* why a->cols + cols - 2 * (b->cols & ~1) ?
651 * what we really want is ceiling((a->cols - (b->cols & ~1)) / (cols - (b->cols & ~1)))
652 * in this case, we strip out paddings on the left/right, and compute how many tiles
653 * we need. It then be interpreted in the above integer division form */
654 int tile_x = ccv_max(1, (a->cols + cols - 2 * (b->cols & ~1)) / (cols - (b->cols & ~1)))({ typeof (1) _a = (1); typeof ((a->cols + cols - 2 * (b->
cols & ~1)) / (cols - (b->cols & ~1))) _b = ((a->
cols + cols - 2 * (b->cols & ~1)) / (cols - (b->cols
& ~1))); (_a > _b) ? _a : _b; });
655 int tile_y = ccv_max(1, (a->rows + rows - 2 * (b->rows & ~1)) / (rows - (b->rows & ~1)))({ typeof (1) _a = (1); typeof ((a->rows + rows - 2 * (b->
rows & ~1)) / (rows - (b->rows & ~1))) _b = ((a->
rows + rows - 2 * (b->rows & ~1)) / (rows - (b->rows
& ~1))); (_a > _b) ? _a : _b; });
656 int brows2 = b->rows / 2;
657 int bcols2 = b->cols / 2;
658#define for_block(_for_type, _cpx_type, _for_set, _for_get) \
659 _for_type scale = 1.0 / (rows * cols); \
660 for (i = 0; i < tile_y; i++) \
661 for (j = 0; j < tile_x; j++) \
662 { \
663 int x, y; \
664 memset(fftw_a, 0, rows * cols_2c * ch * sizeof(_for_type)); \
665 int iy = ccv_min(i * (rows - (b->rows & ~1)), ccv_max(a->rows - rows, 0))({ typeof (i * (rows - (b->rows & ~1))) _a = (i * (rows
- (b->rows & ~1))); typeof (({ typeof (a->rows - rows
) _a = (a->rows - rows); typeof (0) _b = (0); (_a > _b)
? _a : _b; })) _b = (({ typeof (a->rows - rows) _a = (a->
rows - rows); typeof (0) _b = (0); (_a > _b) ? _a : _b; })
); (_a < _b) ? _a : _b; }); \
666 int ix = ccv_min(j * (cols - (b->cols & ~1)), ccv_max(a->cols - cols, 0))({ typeof (j * (cols - (b->cols & ~1))) _a = (j * (cols
- (b->cols & ~1))); typeof (({ typeof (a->cols - cols
) _a = (a->cols - cols); typeof (0) _b = (0); (_a > _b)
? _a : _b; })) _b = (({ typeof (a->cols - cols) _a = (a->
cols - cols); typeof (0) _b = (0); (_a > _b) ? _a : _b; })
); (_a < _b) ? _a : _b; }); \
667 _for_type* fftw_ptr = (_for_type*)fftw_a; \
668 int end_y = ccv_min(rows, a->rows - iy)({ typeof (rows) _a = (rows); typeof (a->rows - iy) _b = (
a->rows - iy); (_a < _b) ? _a : _b; }); \
669 int end_x = ccv_min(cols, a->cols - ix)({ typeof (cols) _a = (cols); typeof (a->cols - ix) _b = (
a->cols - ix); (_a < _b) ? _a : _b; }); \
670 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(a, iy, ix, 0)((((a)->type) & CCV_32S) ? (void*)((a)->data.i32 + (
(iy) * (a)->cols + (ix)) * (((a)->type) & 0xFFF) + (
0)) : ((((a)->type) & CCV_32F) ? (void*)((a)->data.
f32+ ((iy) * (a)->cols + (ix)) * (((a)->type) & 0xFFF
) + (0)) : ((((a)->type) & CCV_64S) ? (void*)((a)->
data.i64+ ((iy) * (a)->cols + (ix)) * (((a)->type) &
0xFFF) + (0)) : ((((a)->type) & CCV_64F) ? (void*)((a
)->data.f64 + ((iy) * (a)->cols + (ix)) * (((a)->type
) & 0xFFF) + (0)) : (void*)((a)->data.u8 + (iy) * (a)->
step + (ix) * (((a)->type) & 0xFFF) + (0)))))); \
671 for (y = 0; y < end_y; y++) \
672 { \
673 for (x = 0; x < end_x * ch; x++) \
674 fftw_ptr[x] = _for_get(m_ptr, x); \
675 fftw_ptr += cols_2c * ch; \
676 m_ptr += a->step; \
677 } \
678 _cpx_type* fftw_ac = (_cpx_type*)fftw_a; \
679 _cpx_type* fftw_bc = (_cpx_type*)fftw_b; \
680 _cpx_type* fftw_dc = (_cpx_type*)fftw_d; \
681 fft_execute_dft_r2c((_for_type*)fftw_a, (_cpx_type*)fftw_ac); \
682 for (x = 0; x < rows * ch * (cols / 2 + 1); x++) \
683 fftw_dc[x] = (fftw_ac[x] * fftw_bc[x]) * scale; \
684 fft_execute_dft_c2r((_cpx_type*)fftw_dc, (_for_type*)fftw_d); \
685 fftw_ptr = (_for_type*)fftw_d + ((1 + (i > 0)) * brows2 * cols_2c + (1 + (j > 0)) * bcols2) * ch; \
686 end_y = ccv_min(d->rows - (iy + (i > 0) * brows2), \({ typeof (d->rows - (iy + (i > 0) * brows2)) _a = (d->
rows - (iy + (i > 0) * brows2)); typeof ((rows - (b->rows
& ~1)) + (i == 0) * brows2) _b = ((rows - (b->rows &
~1)) + (i == 0) * brows2); (_a < _b) ? _a : _b; })
687 (rows - (b->rows & ~1)) + (i == 0) * brows2)({ typeof (d->rows - (iy + (i > 0) * brows2)) _a = (d->
rows - (iy + (i > 0) * brows2)); typeof ((rows - (b->rows
& ~1)) + (i == 0) * brows2) _b = ((rows - (b->rows &
~1)) + (i == 0) * brows2); (_a < _b) ? _a : _b; }); \
688 end_x = ccv_min(d->cols - (ix + (j > 0) * bcols2), \({ typeof (d->cols - (ix + (j > 0) * bcols2)) _a = (d->
cols - (ix + (j > 0) * bcols2)); typeof ((cols - (b->cols
& ~1)) + (j == 0) * bcols2) _b = ((cols - (b->cols &
~1)) + (j == 0) * bcols2); (_a < _b) ? _a : _b; })
689 (cols - (b->cols & ~1)) + (j == 0) * bcols2)({ typeof (d->cols - (ix + (j > 0) * bcols2)) _a = (d->
cols - (ix + (j > 0) * bcols2)); typeof ((cols - (b->cols
& ~1)) + (j == 0) * bcols2) _b = ((cols - (b->cols &
~1)) + (j == 0) * bcols2); (_a < _b) ? _a : _b; }); \
690 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2, ix + (j > 0) * bcols2, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2) * (d)->cols + (ix + (j > 0) *
bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy + (i >
0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2)) * (
((d)->type) & 0xFFF) + (0)) : ((((d)->type) & CCV_64S
) ? (void*)((d)->data.i64+ ((iy + (i > 0) * brows2) * (
d)->cols + (ix + (j > 0) * bcols2)) * (((d)->type) &
0xFFF) + (0)) : ((((d)->type) & CCV_64F) ? (void*)((d
)->data.f64 + ((iy + (i > 0) * brows2) * (d)->cols +
(ix + (j > 0) * bcols2)) * (((d)->type) & 0xFFF) +
(0)) : (void*)((d)->data.u8 + (iy + (i > 0) * brows2) *
(d)->step + (ix + (j > 0) * bcols2) * (((d)->type) &
0xFFF) + (0)))))); \
691 for (y = 0; y < end_y; y++) \
692 { \
693 for (x = 0; x < end_x * ch; x++) \
694 _for_set(m_ptr, x, fftw_ptr[x]); \
695 m_ptr += d->step; \
696 fftw_ptr += cols_2c * ch; \
697 } \
698 int end_tile_y, end_tile_x; \
699 /* handle edge cases: */ \
700 if (i + 1 == tile_y && end_y + iy + (i > 0) * brows2 < d->rows) \
701 { \
702 end_tile_y = ccv_min(brows2, d->rows - (iy + (i > 0) * brows2 + end_y))({ typeof (brows2) _a = (brows2); typeof (d->rows - (iy + (
i > 0) * brows2 + end_y)) _b = (d->rows - (iy + (i >
0) * brows2 + end_y)); (_a < _b) ? _a : _b; }); \
703 fftw_ptr = (_for_type*)fftw_d + (1 + (j > 0)) * bcols2 * ch; \
704 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2 + end_y, ix + (j > 0) * bcols2, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)
->type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy + (
i > 0) * brows2 + end_y) * (d)->cols + (ix + (j > 0)
* bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_64S) ? (void*)((d)->data.i64+ ((iy + (i >
0) * brows2 + end_y) * (d)->cols + (ix + (j > 0) * bcols2
)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->type) &
CCV_64F) ? (void*)((d)->data.f64 + ((iy + (i > 0) * brows2
+ end_y) * (d)->cols + (ix + (j > 0) * bcols2)) * (((d
)->type) & 0xFFF) + (0)) : (void*)((d)->data.u8 + (
iy + (i > 0) * brows2 + end_y) * (d)->step + (ix + (j >
0) * bcols2) * (((d)->type) & 0xFFF) + (0)))))); \
705 for (y = 0; y < end_tile_y; y++) \
706 { \
707 for (x = 0; x < end_x * ch; x++) \
708 _for_set(m_ptr, x, fftw_ptr[x]); \
709 m_ptr += d->step; \
710 fftw_ptr += cols_2c * ch; \
711 } \
712 } \
713 if (j + 1 == tile_x && end_x + ix + (j > 0) * bcols2 < d->cols) \
714 { \
715 end_tile_x = ccv_min(bcols2, d->cols - (ix + (j > 0) * bcols2 + end_x))({ typeof (bcols2) _a = (bcols2); typeof (d->cols - (ix + (
j > 0) * bcols2 + end_x)) _b = (d->cols - (ix + (j >
0) * bcols2 + end_x)); (_a < _b) ? _a : _b; }); \
716 fftw_ptr = (_for_type*)fftw_d + (1 + (i > 0)) * brows2 * cols_2c * ch; \
717 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2, ix + (j > 0) * bcols2 + end_x, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2) * (d)->cols + (ix + (j > 0) *
bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) : (((
(d)->type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy
+ (i > 0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2
+ end_x)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_64S) ? (void*)((d)->data.i64+ ((iy + (i >
0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2 + end_x
)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->type) &
CCV_64F) ? (void*)((d)->data.f64 + ((iy + (i > 0) * brows2
) * (d)->cols + (ix + (j > 0) * bcols2 + end_x)) * (((d
)->type) & 0xFFF) + (0)) : (void*)((d)->data.u8 + (
iy + (i > 0) * brows2) * (d)->step + (ix + (j > 0) *
bcols2 + end_x) * (((d)->type) & 0xFFF) + (0)))))); \
718 for (y = 0; y < end_y; y++) \
719 { \
720 for (x = 0; x < end_tile_x * ch; x++) \
721 _for_set(m_ptr, x, fftw_ptr[x]); \
722 m_ptr += d->step; \
723 fftw_ptr += cols_2c * ch; \
724 } \
725 } \
726 if (i + 1 == tile_y && end_y + iy + (i > 0) * brows2 < d->rows && \
727 j + 1 == tile_x && end_x + ix + (j > 0) * bcols2 < d->cols) \
728 { \
729 fftw_ptr = (_for_type*)fftw_d; \
730 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2 + end_y, ix + (j > 0) * bcols2 + end_x, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_32F) ? (void*)((d)->data.f32+ (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_64S) ? (void*)((d)->data.i64+ (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_64F) ? (void*)((d)->data.f64 +
((iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (
j > 0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) +
(0)) : (void*)((d)->data.u8 + (iy + (i > 0) * brows2 +
end_y) * (d)->step + (ix + (j > 0) * bcols2 + end_x) *
(((d)->type) & 0xFFF) + (0)))))); \
731 for (y = 0; y < end_tile_y; y++) \
732 { \
733 for (x = 0; x < end_tile_x * ch; x++) \
734 _for_set(m_ptr, x, fftw_ptr[x]); \
735 m_ptr += d->step; \
736 fftw_ptr += cols_2c * ch; \
737 } \
738 } \
739 }
740 if (fft_type == CCV_32F)
741 {
742#define fft_execute_dft_r2c(r, c) fftwf_execute_dft_r2c(pf, r, c)
743#define fft_execute_dft_c2r(c, r) fftwf_execute_dft_c2r(pinvf, c, r)
744 ccv_matrix_setter(d->type, ccv_matrix_getter, a->type, for_block, float, fftwf_complex){ switch (((d->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(float
, fftwf_complex, _ccv_set_32s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(float, fftwf_complex, _ccv_set_32s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(float
, fftwf_complex, _ccv_set_32s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(float, fftwf_complex, _ccv_set_32s_value
, _ccv_get_64f_value); break; } default: { for_block(float, fftwf_complex
, _ccv_set_32s_value, _ccv_get_8u_value); } } }; break; } case
CCV_32F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(float, fftwf_complex, _ccv_set_32f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, fftwf_complex, _ccv_set_32f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(float
, fftwf_complex, _ccv_set_32f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(float, fftwf_complex, _ccv_set_32f_value
, _ccv_get_64f_value); break; } default: { for_block(float, fftwf_complex
, _ccv_set_32f_value, _ccv_get_8u_value); } } }; break; } case
CCV_64S: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(float, fftwf_complex, _ccv_set_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, fftwf_complex, _ccv_set_64s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(float
, fftwf_complex, _ccv_set_64s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(float, fftwf_complex, _ccv_set_64s_value
, _ccv_get_64f_value); break; } default: { for_block(float, fftwf_complex
, _ccv_set_64s_value, _ccv_get_8u_value); } } }; break; } case
CCV_64F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(float, fftwf_complex, _ccv_set_64f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, fftwf_complex, _ccv_set_64f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(float
, fftwf_complex, _ccv_set_64f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(float, fftwf_complex, _ccv_set_64f_value
, _ccv_get_64f_value); break; } default: { for_block(float, fftwf_complex
, _ccv_set_64f_value, _ccv_get_8u_value); } } }; break; } default
: { { switch (((a->type) & 0xFF000)) { case CCV_32S: {
for_block(float, fftwf_complex, _ccv_set_8u_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, fftwf_complex, _ccv_set_8u_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(float
, fftwf_complex, _ccv_set_8u_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(float, fftwf_complex, _ccv_set_8u_value
, _ccv_get_64f_value); break; } default: { for_block(float, fftwf_complex
, _ccv_set_8u_value, _ccv_get_8u_value); } } }; } } };
745#undef fft_execute_dft_r2c
746#undef fft_execute_dft_c2r
747 pthread_mutex_lock(&fftw_plan_mutex);
748 fftwf_destroy_plan(pf);
749 fftwf_destroy_plan(pinvf);
750 pthread_mutex_unlock(&fftw_plan_mutex);
751 fftwf_free(fftw_a);
752 fftwf_free(fftw_b);
753 fftwf_free(fftw_d);
754 } else {
755#define fft_execute_dft_r2c(r, c) fftw_execute_dft_r2c(p, r, c)
756#define fft_execute_dft_c2r(c, r) fftw_execute_dft_c2r(pinv, c, r)
757 ccv_matrix_setter(d->type, ccv_matrix_getter, a->type, for_block, double, fftw_complex){ switch (((d->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(double
, fftw_complex, _ccv_set_32s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(double, fftw_complex, _ccv_set_32s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(double
, fftw_complex, _ccv_set_32s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, fftw_complex, _ccv_set_32s_value
, _ccv_get_64f_value); break; } default: { for_block(double, fftw_complex
, _ccv_set_32s_value, _ccv_get_8u_value); } } }; break; } case
CCV_32F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(double, fftw_complex, _ccv_set_32f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(double, fftw_complex, _ccv_set_32f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(double
, fftw_complex, _ccv_set_32f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, fftw_complex, _ccv_set_32f_value
, _ccv_get_64f_value); break; } default: { for_block(double, fftw_complex
, _ccv_set_32f_value, _ccv_get_8u_value); } } }; break; } case
CCV_64S: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(double, fftw_complex, _ccv_set_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(double, fftw_complex, _ccv_set_64s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(double
, fftw_complex, _ccv_set_64s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, fftw_complex, _ccv_set_64s_value
, _ccv_get_64f_value); break; } default: { for_block(double, fftw_complex
, _ccv_set_64s_value, _ccv_get_8u_value); } } }; break; } case
CCV_64F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(double, fftw_complex, _ccv_set_64f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(double, fftw_complex, _ccv_set_64f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(double
, fftw_complex, _ccv_set_64f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, fftw_complex, _ccv_set_64f_value
, _ccv_get_64f_value); break; } default: { for_block(double, fftw_complex
, _ccv_set_64f_value, _ccv_get_8u_value); } } }; break; } default
: { { switch (((a->type) & 0xFF000)) { case CCV_32S: {
for_block(double, fftw_complex, _ccv_set_8u_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(double, fftw_complex, _ccv_set_8u_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(double
, fftw_complex, _ccv_set_8u_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, fftw_complex, _ccv_set_8u_value
, _ccv_get_64f_value); break; } default: { for_block(double, fftw_complex
, _ccv_set_8u_value, _ccv_get_8u_value); } } }; } } };
758#undef fft_execute_dft_r2c
759#undef fft_execute_dft_c2r
760 pthread_mutex_lock(&fftw_plan_mutex);
761 fftw_destroy_plan(p);
762 fftw_destroy_plan(pinv);
763 pthread_mutex_unlock(&fftw_plan_mutex);
764 fftw_free(fftw_a);
765 fftw_free(fftw_b);
766 fftw_free(fftw_d);
767 }
768#undef for_block
769}
770#else
771static void _ccv_filter_kissfft(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t* d, int padding_pattern)
772{
773 int ch = CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF);
774 int fft_type = (CCV_GET_DATA_TYPE(d->type)((d->type) & 0xFF000) == CCV_8U || CCV_GET_DATA_TYPE(d->type)((d->type) & 0xFF000) == CCV_32F) ? CCV_32F : CCV_64F;
775 int rows = ((ccv_min(a->rows + b->rows - 1, kiss_fftr_next_fast_size_real(b->rows * 3))({ typeof (a->rows + b->rows - 1) _a = (a->rows + b->
rows - 1); typeof (kiss_fftr_next_fast_size_real(b->rows *
3)) _b = (kiss_fftr_next_fast_size_real(b->rows * 3)); (_a
< _b) ? _a : _b; }) + 1) >> 1) << 1;
776 int cols = ((ccv_min(a->cols + b->cols - 1, kiss_fftr_next_fast_size_real(b->cols * 3))({ typeof (a->cols + b->cols - 1) _a = (a->cols + b->
cols - 1); typeof (kiss_fftr_next_fast_size_real(b->cols *
3)) _b = (kiss_fftr_next_fast_size_real(b->cols * 3)); (_a
< _b) ? _a : _b; }) + 1) >> 1) << 1;
777 int ndim[] = {rows, cols};
778 void* kiss_a;
779 void* kiss_b;
780 void* kiss_d;
781 void* kiss_ac;
782 void* kiss_bc;
783 void* kiss_dc;
784 kiss_fftndr_cfg p;
785 kiss_fftndr_cfg pinv;
786 kissf_fftndr_cfg pf;
787 kissf_fftndr_cfg pinvf;
788 if (fft_type == CCV_32F)
789 {
790 pf = kissf_fftndr_alloc(ndim, 2, 0, 0, 0);
791 pinvf = kissf_fftndr_alloc(ndim, 2, 1, 0, 0);
792 kiss_a = ccmallocmalloc(rows * cols * ch * sizeof(kissf_fft_scalar));
793 kiss_b = ccmallocmalloc(rows * cols * ch * sizeof(kissf_fft_scalar));
794 kiss_d = ccmallocmalloc(rows * cols * ch * sizeof(kissf_fft_scalar));
795 kiss_ac = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kissf_fft_cpx));
796 kiss_bc = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kissf_fft_cpx));
797 kiss_dc = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kissf_fft_cpx));
798 memset(kiss_b, 0, rows * cols * ch * sizeof(kissf_fft_scalar));
799 } else {
800 p = kiss_fftndr_alloc(ndim, 2, 0, 0, 0);
801 pinv = kiss_fftndr_alloc(ndim, 2, 1, 0, 0);
802 kiss_a = ccmallocmalloc(rows * cols * ch * sizeof(kiss_fft_scalar));
803 kiss_b = ccmallocmalloc(rows * cols * ch * sizeof(kiss_fft_scalar));
804 kiss_d = ccmallocmalloc(rows * cols * ch * sizeof(kiss_fft_scalar));
805 kiss_ac = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kiss_fft_cpx));
806 kiss_bc = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kiss_fft_cpx));
807 kiss_dc = ccmallocmalloc(rows * (cols / 2 + 1) * ch * sizeof(kiss_fft_cpx));
808 memset(kiss_b, 0, rows * cols * ch * sizeof(kiss_fft_scalar));
809 }
810 int nch = rows * cols, nchc = rows * (cols / 2 + 1);
811 int i, j, k;
812 unsigned char* m_ptr = b->data.u8;
813 // to flip matrix b is crucial, this problem only shows when I changed to a more sophisticated test case
814#define for_block(_for_type, _for_get) \
815 _for_type* kiss_ptr = (_for_type*)kiss_b + (b->rows - 1) * cols; \
816 for (i = 0; i < b->rows; i++) \
817 { \
818 for (j = 0; j < b->cols; j++) \
819 for (k = 0; k < ch; k++) \
820 kiss_ptr[k * nch + b->cols - 1 - j] = _for_get(m_ptr, j * ch + k); \
821 kiss_ptr -= cols; \
822 m_ptr += b->step; \
823 }
824 ccv_matrix_typeof(fft_type, ccv_matrix_getter, b->type, for_block){ switch (((fft_type) & 0xFF000)) { case CCV_32S: { { switch
(((b->type) & 0xFF000)) { case CCV_32S: { for_block(int
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(int
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(int
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(int
, _ccv_get_64f_value); break; } default: { for_block(int, _ccv_get_8u_value
); } } }; break; } case CCV_32F: { { switch (((b->type) &
0xFF000)) { case CCV_32S: { for_block(float, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(float, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(float, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(float, _ccv_get_64f_value
); break; } default: { for_block(float, _ccv_get_8u_value); }
} }; break; } case CCV_64S: { { switch (((b->type) & 0xFF000
)) { case CCV_32S: { for_block(int64_t, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(int64_t, _ccv_get_32f_value); break
; } case CCV_64S: { for_block(int64_t, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(int64_t, _ccv_get_64f_value); break
; } default: { for_block(int64_t, _ccv_get_8u_value); } } }; break
; } case CCV_64F: { { switch (((b->type) & 0xFF000)) {
case CCV_32S: { for_block(double, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(double, _ccv_get_32f_value); break
; } case CCV_64S: { for_block(double, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(double, _ccv_get_64f_value); break
; } default: { for_block(double, _ccv_get_8u_value); } } }; break
; } default: { { switch (((b->type) & 0xFF000)) { case
CCV_32S: { for_block(unsigned char, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(unsigned char, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(unsigned char, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(unsigned char, _ccv_get_64f_value
); break; } default: { for_block(unsigned char, _ccv_get_8u_value
); } } }; } } };
825#undef for_block
826 if (fft_type == CCV_32F)
827 for (k = 0; k < ch; k++)
828 kissf_fftndr(pf, (kissf_fft_scalar*)kiss_b + nch * k, (kissf_fft_cpx*)kiss_bc + nchc * k);
829 else
830 for (k = 0; k < ch; k++)
831 kiss_fftndr(p, (kiss_fft_scalar*)kiss_b + nch * k, (kiss_fft_cpx*)kiss_bc + nchc * k);
832 /* why a->cols + cols - 2 * (b->cols & ~1) ?
833 * what we really want is ceiling((a->cols - (b->cols & ~1)) / (cols - (b->cols & ~1)))
834 * in this case, we strip out paddings on the left/right, and compute how many tiles
835 * we need. It then be interpreted in the above integer division form */
836 int tile_x = ccv_max(1, (a->cols + cols - 2 * (b->cols & ~1)) / (cols - (b->cols & ~1)))({ typeof (1) _a = (1); typeof ((a->cols + cols - 2 * (b->
cols & ~1)) / (cols - (b->cols & ~1))) _b = ((a->
cols + cols - 2 * (b->cols & ~1)) / (cols - (b->cols
& ~1))); (_a > _b) ? _a : _b; });
837 int tile_y = ccv_max(1, (a->rows + rows - 2 * (b->rows & ~1)) / (rows - (b->rows & ~1)))({ typeof (1) _a = (1); typeof ((a->rows + rows - 2 * (b->
rows & ~1)) / (rows - (b->rows & ~1))) _b = ((a->
rows + rows - 2 * (b->rows & ~1)) / (rows - (b->rows
& ~1))); (_a > _b) ? _a : _b; });
838 int brows2 = b->rows / 2;
839 int bcols2 = b->cols / 2;
840#define for_block(_for_type, _cpx_type, _for_set, _for_get) \
841 _for_type scale = 1.0 / (rows * cols); \
842 for (i = 0; i < tile_y; i++) \
843 for (j = 0; j < tile_x; j++) \
844 { \
845 int x, y; \
846 memset(kiss_a, 0, rows * cols * ch * sizeof(_for_type)); \
847 int iy = ccv_min(i * (rows - (b->rows & ~1)), ccv_max(a->rows - rows, 0))({ typeof (i * (rows - (b->rows & ~1))) _a = (i * (rows
- (b->rows & ~1))); typeof (({ typeof (a->rows - rows
) _a = (a->rows - rows); typeof (0) _b = (0); (_a > _b)
? _a : _b; })) _b = (({ typeof (a->rows - rows) _a = (a->
rows - rows); typeof (0) _b = (0); (_a > _b) ? _a : _b; })
); (_a < _b) ? _a : _b; }); \
848 int ix = ccv_min(j * (cols - (b->cols & ~1)), ccv_max(a->cols - cols, 0))({ typeof (j * (cols - (b->cols & ~1))) _a = (j * (cols
- (b->cols & ~1))); typeof (({ typeof (a->cols - cols
) _a = (a->cols - cols); typeof (0) _b = (0); (_a > _b)
? _a : _b; })) _b = (({ typeof (a->cols - cols) _a = (a->
cols - cols); typeof (0) _b = (0); (_a > _b) ? _a : _b; })
); (_a < _b) ? _a : _b; }); \
849 _for_type* kiss_ptr = (_for_type*)kiss_a; \
850 int end_y = ccv_min(rows, a->rows - iy)({ typeof (rows) _a = (rows); typeof (a->rows - iy) _b = (
a->rows - iy); (_a < _b) ? _a : _b; }); \
851 int end_x = ccv_min(cols, a->cols - ix)({ typeof (cols) _a = (cols); typeof (a->cols - ix) _b = (
a->cols - ix); (_a < _b) ? _a : _b; }); \
852 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(a, iy, ix, 0)((((a)->type) & CCV_32S) ? (void*)((a)->data.i32 + (
(iy) * (a)->cols + (ix)) * (((a)->type) & 0xFFF) + (
0)) : ((((a)->type) & CCV_32F) ? (void*)((a)->data.
f32+ ((iy) * (a)->cols + (ix)) * (((a)->type) & 0xFFF
) + (0)) : ((((a)->type) & CCV_64S) ? (void*)((a)->
data.i64+ ((iy) * (a)->cols + (ix)) * (((a)->type) &
0xFFF) + (0)) : ((((a)->type) & CCV_64F) ? (void*)((a
)->data.f64 + ((iy) * (a)->cols + (ix)) * (((a)->type
) & 0xFFF) + (0)) : (void*)((a)->data.u8 + (iy) * (a)->
step + (ix) * (((a)->type) & 0xFFF) + (0)))))); \
853 for (y = 0; y < end_y; y++) \
854 { \
855 for (x = 0; x < end_x; x++) \
856 for (k = 0; k < ch; k++) \
857 kiss_ptr[k * nch + x] = _for_get(m_ptr, x * ch + k); \
858 kiss_ptr += cols; \
859 m_ptr += a->step; \
860 } \
861 for (k = 0; k < ch; k++) \
862 fft_ndr((_for_type*)kiss_a + nch * k, (_cpx_type*)kiss_ac + nchc * k); \
863 _cpx_type* fft_ac = (_cpx_type*)kiss_ac; \
864 _cpx_type* fft_bc = (_cpx_type*)kiss_bc; \
865 _cpx_type* fft_dc = (_cpx_type*)kiss_dc; \
866 for (x = 0; x < rows * ch * (cols / 2 + 1); x++) \
867 { \
868 fft_dc[x].r = (fft_ac[x].r * fft_bc[x].r - fft_ac[x].i * fft_bc[x].i) * scale; \
869 fft_dc[x].i = (fft_ac[x].i * fft_bc[x].r + fft_ac[x].r * fft_bc[x].i) * scale; \
870 } \
871 for (k = 0; k < ch; k++) \
872 fft_ndri((_cpx_type*)kiss_dc + nchc * k, (_for_type*)kiss_d + nch * k); \
873 kiss_ptr = (_for_type*)kiss_d + (1 + (i > 0)) * brows2 * cols + (1 + (j > 0)) * bcols2; \
874 end_y = ccv_min(d->rows - (iy + (i > 0) * brows2), \({ typeof (d->rows - (iy + (i > 0) * brows2)) _a = (d->
rows - (iy + (i > 0) * brows2)); typeof ((rows - (b->rows
& ~1)) + (i == 0) * brows2) _b = ((rows - (b->rows &
~1)) + (i == 0) * brows2); (_a < _b) ? _a : _b; })
875 (rows - (b->rows & ~1)) + (i == 0) * brows2)({ typeof (d->rows - (iy + (i > 0) * brows2)) _a = (d->
rows - (iy + (i > 0) * brows2)); typeof ((rows - (b->rows
& ~1)) + (i == 0) * brows2) _b = ((rows - (b->rows &
~1)) + (i == 0) * brows2); (_a < _b) ? _a : _b; }); \
876 end_x = ccv_min(d->cols - (ix + (j > 0) * bcols2), \({ typeof (d->cols - (ix + (j > 0) * bcols2)) _a = (d->
cols - (ix + (j > 0) * bcols2)); typeof ((cols - (b->cols
& ~1)) + (j == 0) * bcols2) _b = ((cols - (b->cols &
~1)) + (j == 0) * bcols2); (_a < _b) ? _a : _b; })
877 (cols - (b->cols & ~1)) + (j == 0) * bcols2)({ typeof (d->cols - (ix + (j > 0) * bcols2)) _a = (d->
cols - (ix + (j > 0) * bcols2)); typeof ((cols - (b->cols
& ~1)) + (j == 0) * bcols2) _b = ((cols - (b->cols &
~1)) + (j == 0) * bcols2); (_a < _b) ? _a : _b; }); \
878 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2, ix + (j > 0) * bcols2, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2) * (d)->cols + (ix + (j > 0) *
bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy + (i >
0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2)) * (
((d)->type) & 0xFFF) + (0)) : ((((d)->type) & CCV_64S
) ? (void*)((d)->data.i64+ ((iy + (i > 0) * brows2) * (
d)->cols + (ix + (j > 0) * bcols2)) * (((d)->type) &
0xFFF) + (0)) : ((((d)->type) & CCV_64F) ? (void*)((d
)->data.f64 + ((iy + (i > 0) * brows2) * (d)->cols +
(ix + (j > 0) * bcols2)) * (((d)->type) & 0xFFF) +
(0)) : (void*)((d)->data.u8 + (iy + (i > 0) * brows2) *
(d)->step + (ix + (j > 0) * bcols2) * (((d)->type) &
0xFFF) + (0)))))); \
879 for (y = 0; y < end_y; y++) \
880 { \
881 for (x = 0; x < end_x; x++) \
882 for (k = 0; k < ch; k++) \
883 _for_set(m_ptr, x * ch + k, kiss_ptr[k * nch + x]); \
884 m_ptr += d->step; \
885 kiss_ptr += cols; \
886 } \
887 int end_tile_y, end_tile_x; \
888 /* handle edge cases: */ \
889 if (i + 1 == tile_y && end_y + iy + (i > 0) * brows2 < d->rows) \
890 { \
891 end_tile_y = ccv_min(brows2, d->rows - (iy + (i > 0) * brows2 + end_y))({ typeof (brows2) _a = (brows2); typeof (d->rows - (iy + (
i > 0) * brows2 + end_y)) _b = (d->rows - (iy + (i >
0) * brows2 + end_y)); (_a < _b) ? _a : _b; }); \
892 kiss_ptr = (_for_type*)kiss_d + (1 + (j > 0)) * bcols2; \
893 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2 + end_y, ix + (j > 0) * bcols2, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)
->type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy + (
i > 0) * brows2 + end_y) * (d)->cols + (ix + (j > 0)
* bcols2)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_64S) ? (void*)((d)->data.i64+ ((iy + (i >
0) * brows2 + end_y) * (d)->cols + (ix + (j > 0) * bcols2
)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->type) &
CCV_64F) ? (void*)((d)->data.f64 + ((iy + (i > 0) * brows2
+ end_y) * (d)->cols + (ix + (j > 0) * bcols2)) * (((d
)->type) & 0xFFF) + (0)) : (void*)((d)->data.u8 + (
iy + (i > 0) * brows2 + end_y) * (d)->step + (ix + (j >
0) * bcols2) * (((d)->type) & 0xFFF) + (0)))))); \
894 for (y = 0; y < end_tile_y; y++) \
895 { \
896 for (x = 0; x < end_x; x++) \
897 for (k = 0; k < ch; k++) \
898 _for_set(m_ptr, x * ch + k, kiss_ptr[k * nch + x]); \
899 m_ptr += d->step; \
900 kiss_ptr += cols; \
901 } \
902 } \
903 if (j + 1 == tile_x && end_x + ix + (j > 0) * bcols2 < d->cols) \
904 { \
905 end_tile_x = ccv_min(bcols2, d->cols - (ix + (j > 0) * bcols2 + end_x))({ typeof (bcols2) _a = (bcols2); typeof (d->cols - (ix + (
j > 0) * bcols2 + end_x)) _b = (d->cols - (ix + (j >
0) * bcols2 + end_x)); (_a < _b) ? _a : _b; }); \
906 kiss_ptr = (_for_type*)kiss_d + (1 + (i > 0)) * brows2 * cols; \
907 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2, ix + (j > 0) * bcols2 + end_x, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2) * (d)->cols + (ix + (j > 0) *
bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) : (((
(d)->type) & CCV_32F) ? (void*)((d)->data.f32+ ((iy
+ (i > 0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2
+ end_x)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->
type) & CCV_64S) ? (void*)((d)->data.i64+ ((iy + (i >
0) * brows2) * (d)->cols + (ix + (j > 0) * bcols2 + end_x
)) * (((d)->type) & 0xFFF) + (0)) : ((((d)->type) &
CCV_64F) ? (void*)((d)->data.f64 + ((iy + (i > 0) * brows2
) * (d)->cols + (ix + (j > 0) * bcols2 + end_x)) * (((d
)->type) & 0xFFF) + (0)) : (void*)((d)->data.u8 + (
iy + (i > 0) * brows2) * (d)->step + (ix + (j > 0) *
bcols2 + end_x) * (((d)->type) & 0xFFF) + (0)))))); \
908 for (y = 0; y < end_y; y++) \
909 { \
910 for (x = 0; x < end_tile_x; x++) \
911 for (k = 0; k < ch; k++) \
912 _for_set(m_ptr, x * ch + k, kiss_ptr[k * nch + x]); \
913 m_ptr += d->step; \
914 kiss_ptr += cols; \
915 } \
916 } \
917 if (i + 1 == tile_y && end_y + iy + (i > 0) * brows2 < d->rows && \
918 j + 1 == tile_x && end_x + ix + (j > 0) * bcols2 < d->cols) \
919 { \
920 kiss_ptr = (_for_type*)kiss_d; \
921 m_ptr = (unsigned char*)ccv_get_dense_matrix_cell(d, iy + (i > 0) * brows2 + end_y, ix + (j > 0) * bcols2 + end_x, 0)((((d)->type) & CCV_32S) ? (void*)((d)->data.i32 + (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_32F) ? (void*)((d)->data.f32+ (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_64S) ? (void*)((d)->data.i64+ (
(iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (j >
0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) + (0)) :
((((d)->type) & CCV_64F) ? (void*)((d)->data.f64 +
((iy + (i > 0) * brows2 + end_y) * (d)->cols + (ix + (
j > 0) * bcols2 + end_x)) * (((d)->type) & 0xFFF) +
(0)) : (void*)((d)->data.u8 + (iy + (i > 0) * brows2 +
end_y) * (d)->step + (ix + (j > 0) * bcols2 + end_x) *
(((d)->type) & 0xFFF) + (0)))))); \
922 for (y = 0; y < end_tile_y; y++) \
923 { \
924 for (x = 0; x < end_tile_x; x++) \
925 for (k = 0; k < ch; k++) \
926 _for_set(m_ptr, x * ch + k, kiss_ptr[k * nch + x]); \
927 m_ptr += d->step; \
928 kiss_ptr += cols; \
929 } \
930 } \
931 }
932 if (fft_type == CCV_32F)
933 {
934#define fft_ndr(r, c) kissf_fftndr(pf, r, c)
935#define fft_ndri(c, r) kissf_fftndri(pinvf, c, r)
936 ccv_matrix_setter(d->type, ccv_matrix_getter, a->type, for_block, kissf_fft_scalar, kissf_fft_cpx){ switch (((d->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_32s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kissf_fft_scalar, kissf_fft_cpx
, _ccv_set_32s_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_32s_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_32s_value, _ccv_get_64f_value); break
; } default: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_32s_value
, _ccv_get_8u_value); } } }; break; } case CCV_32F: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_32f_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kissf_fft_scalar, kissf_fft_cpx
, _ccv_set_32f_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_32f_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_32f_value, _ccv_get_64f_value); break
; } default: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_32f_value
, _ccv_get_8u_value); } } }; break; } case CCV_64S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_64s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kissf_fft_scalar, kissf_fft_cpx
, _ccv_set_64s_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_64s_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_64s_value, _ccv_get_64f_value); break
; } default: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_64s_value
, _ccv_get_8u_value); } } }; break; } case CCV_64F: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_64f_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kissf_fft_scalar, kissf_fft_cpx
, _ccv_set_64f_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_64f_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_64f_value, _ccv_get_64f_value); break
; } default: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_64f_value
, _ccv_get_8u_value); } } }; break; } default: { { switch (((
a->type) & 0xFF000)) { case CCV_32S: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_8u_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kissf_fft_scalar, kissf_fft_cpx
, _ccv_set_8u_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_8u_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kissf_fft_scalar
, kissf_fft_cpx, _ccv_set_8u_value, _ccv_get_64f_value); break
; } default: { for_block(kissf_fft_scalar, kissf_fft_cpx, _ccv_set_8u_value
, _ccv_get_8u_value); } } }; } } };
937#undef fft_ndr
938#undef fft_ndri
939 kissf_fft_free(pf);
940 kissf_fft_free(pinvf);
941 } else {
942#define fft_ndr(r, c) kiss_fftndr(p, r, c)
943#define fft_ndri(c, r) kiss_fftndri(pinv, c, r)
944 ccv_matrix_setter(d->type, ccv_matrix_getter, a->type, for_block, kiss_fft_scalar, kiss_fft_cpx){ switch (((d->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_32s_value, _ccv_get_32s_value); break
; } case CCV_32F: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_32s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_32s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_32s_value
, _ccv_get_64f_value); break; } default: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_32s_value, _ccv_get_8u_value); } } }
; break; } case CCV_32F: { { switch (((a->type) & 0xFF000
)) { case CCV_32S: { for_block(kiss_fft_scalar, kiss_fft_cpx,
_ccv_set_32f_value, _ccv_get_32s_value); break; } case CCV_32F
: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_32f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_32f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_32f_value
, _ccv_get_64f_value); break; } default: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_32f_value, _ccv_get_8u_value); } } }
; break; } case CCV_64S: { { switch (((a->type) & 0xFF000
)) { case CCV_32S: { for_block(kiss_fft_scalar, kiss_fft_cpx,
_ccv_set_64s_value, _ccv_get_32s_value); break; } case CCV_32F
: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_64s_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_64s_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_64s_value
, _ccv_get_64f_value); break; } default: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_64s_value, _ccv_get_8u_value); } } }
; break; } case CCV_64F: { { switch (((a->type) & 0xFF000
)) { case CCV_32S: { for_block(kiss_fft_scalar, kiss_fft_cpx,
_ccv_set_64f_value, _ccv_get_32s_value); break; } case CCV_32F
: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_64f_value
, _ccv_get_32f_value); break; } case CCV_64S: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_64f_value, _ccv_get_64s_value); break
; } case CCV_64F: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_64f_value
, _ccv_get_64f_value); break; } default: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_64f_value, _ccv_get_8u_value); } } }
; break; } default: { { switch (((a->type) & 0xFF000))
{ case CCV_32S: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_8u_value
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_8u_value, _ccv_get_32f_value); break
; } case CCV_64S: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_8u_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(kiss_fft_scalar
, kiss_fft_cpx, _ccv_set_8u_value, _ccv_get_64f_value); break
; } default: { for_block(kiss_fft_scalar, kiss_fft_cpx, _ccv_set_8u_value
, _ccv_get_8u_value); } } }; } } };
945#undef fft_ndr
946#undef fft_ndri
947 kiss_fft_free(p);
948 kiss_fft_free(pinv);
949 }
950#undef for_block
951 ccfreefree(kiss_dc);
952 ccfreefree(kiss_bc);
953 ccfreefree(kiss_ac);
954 ccfreefree(kiss_d);
955 ccfreefree(kiss_b);
956 ccfreefree(kiss_a);
957}
958#endif
959
960static void _ccv_filter_direct_8u(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t* d, int padding_pattern)
961{
962 int i, j, y, x, k;
963 int nz = b->rows * b->cols;
964 int* coeff = (int*)ccmallocmalloc(nz * sizeof(int));
965 int* cx = (int*)ccmallocmalloc(nz * sizeof(int));
966 int* cy = (int*)ccmallocmalloc(nz * sizeof(int));
967 int scale = 1 << 14;
968 nz = 0;
969 for (i = 0; i < b->rows; i++)
970 for (j = 0; j < b->cols; j++)
971 {
972 coeff[nz] = (int)(ccv_get_dense_matrix_cell_value(b, i, j, 0)((((b)->type) & CCV_32S) ? (b)->data.i32[((i) * (b)
->cols + (j)) * (((b)->type) & 0xFFF) + (0)] : ((((
b)->type) & CCV_32F) ? (b)->data.f32[((i) * (b)->
cols + (j)) * (((b)->type) & 0xFFF) + (0)] : ((((b)->
type) & CCV_64S) ? (b)->data.i64[((i) * (b)->cols +
(j)) * (((b)->type) & 0xFFF) + (0)] : ((((b)->type
) & CCV_64F) ? (b)->data.f64[((i) * (b)->cols + (j)
) * (((b)->type) & 0xFFF) + (0)] : (b)->data.u8[(i)
* (b)->step + (j) * (((b)->type) & 0xFFF) + (0)]))
)) * scale + 0.5);
973 if (coeff[nz] == 0)
974 continue;
975 cy[nz] = i;
976 cx[nz] = j;
977 nz++;
978 }
979 ccv_dense_matrix_t* pa = ccv_dense_matrix_new(a->rows + b->rows / 2 * 2, a->cols + b->cols / 2 * 2, CCV_8U | CCV_C1, 0, 0);
980 /* the padding pattern is different from FFT: |aa{BORDER}|abcd|{BORDER}dd| */
981 for (i = 0; i < pa->rows; i++)
982 for (j = 0; j < pa->cols; j++)
983 pa->data.u8[i * pa->step + j] = a->data.u8[ccv_clamp(i - b->rows / 2, 0, a->rows - 1)({ typeof (0) _a = (0); typeof (a->rows - 1) _b = (a->rows
- 1); typeof (i - b->rows / 2) _x = (i - b->rows / 2);
(_x < _a) ? _a : ((_x > _b) ? _b : _x); }) * a->step + ccv_clamp(j - b->cols / 2, 0, a->cols - 1)({ typeof (0) _a = (0); typeof (a->cols - 1) _b = (a->cols
- 1); typeof (j - b->cols / 2) _x = (j - b->cols / 2);
(_x < _a) ? _a : ((_x > _b) ? _b : _x); })];
984 unsigned char* m_ptr = d->data.u8;
985 unsigned char* a_ptr = pa->data.u8;
986 /* 0.5 denote the overhead for indexing x and y */
987 if (nz < b->rows * b->cols * 0.75)
988 {
989 for (i = 0; i < d->rows; i++)
990 {
991 for (j = 0; j < d->cols; j++)
992 {
993 int z = 0;
994 for (k = 0; k < nz; k++)
995 z += a_ptr[cy[k] * pa->step + j + cx[k]] * coeff[k];
996 m_ptr[j] = ccv_clamp(z >> 14, 0, 255)({ typeof (0) _a = (0); typeof (255) _b = (255); typeof (z >>
14) _x = (z >> 14); (_x < _a) ? _a : ((_x > _b) ?
_b : _x); });
997 }
998 m_ptr += d->step;
999 a_ptr += pa->step;
1000 }
1001 } else {
1002 k = 0;
1003 for (i = 0; i < b->rows; i++)
1004 for (j = 0; j < b->cols; j++)
1005 {
1006 coeff[k] = (int)(ccv_get_dense_matrix_cell_value(b, i, j, 0)((((b)->type) & CCV_32S) ? (b)->data.i32[((i) * (b)
->cols + (j)) * (((b)->type) & 0xFFF) + (0)] : ((((
b)->type) & CCV_32F) ? (b)->data.f32[((i) * (b)->
cols + (j)) * (((b)->type) & 0xFFF) + (0)] : ((((b)->
type) & CCV_64S) ? (b)->data.i64[((i) * (b)->cols +
(j)) * (((b)->type) & 0xFFF) + (0)] : ((((b)->type
) & CCV_64F) ? (b)->data.f64[((i) * (b)->cols + (j)
) * (((b)->type) & 0xFFF) + (0)] : (b)->data.u8[(i)
* (b)->step + (j) * (((b)->type) & 0xFFF) + (0)]))
)) * scale + 0.5);
1007 k++;
1008 }
1009 for (i = 0; i < d->rows; i++)
1010 {
1011 for (j = 0; j < d->cols; j++)
1012 {
1013 int* c_ptr = coeff;
1014 int z = 0;
1015 for (y = 0; y < b->rows; y++)
1016 {
1017 int iyx = y * pa->step;
1018 for (x = 0; x < b->cols; x++)
1019 {
1020 z += a_ptr[iyx + j + x] * c_ptr[0];
1021 c_ptr++;
1022 }
1023 }
1024 m_ptr[j] = ccv_clamp(z >> 14, 0, 255)({ typeof (0) _a = (0); typeof (255) _b = (255); typeof (z >>
14) _x = (z >> 14); (_x < _a) ? _a : ((_x > _b) ?
_b : _x); });
1025 }
1026 m_ptr += d->step;
1027 a_ptr += pa->step;
1028 }
1029 }
1030 ccv_matrix_free(pa);
1031 ccfreefree(coeff);
1032 ccfreefree(cx);
1033 ccfreefree(cy);
1034}
1035
1036void ccv_filter(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t** d, int type, int padding_pattern)
1037{
1038 ccv_declare_derived_signature(sig, a->sig != 0 && b->sig != 0, ccv_sign_with_literal("ccv_filter"), a->sig, b->sig, CCV_EOF_SIGN)char _ccv_identifier_1038[] = ("ccv_filter"); size_t _ccv_string_size_1038
= sizeof(_ccv_identifier_1038);; uint64_t sig = (a->sig !=
0 && b->sig != 0) ? ccv_cache_generate_signature(
_ccv_identifier_1038, _ccv_string_size_1038, a->sig, b->
sig, ((uint64_t)0)) : 0;;
1039 type = (type == 0) ? CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) | CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF) : CCV_GET_DATA_TYPE(type)((type) & 0xFF000) | CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF);
1040 ccv_dense_matrix_t* dd = *d = ccv_dense_matrix_renew(*d, a->rows, a->cols, CCV_ALL_DATA_TYPE(CCV_8U | CCV_32S | CCV_32F | CCV_64S | CCV_64F | CCV_16F | CCV_QX
) | CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF), type, sig);
1041 ccv_object_return_if_cached(, dd){ if ((!(dd) || (((int*)(dd))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE))) { (void)((dd) && (((int*)(dd))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));;
; return ; } };
1042
1043 /* 15 is the constant to indicate the high cost of FFT (even with O(nlog(m)) for
1044 * integer image.
1045 * NOTE: FFT has time complexity of O(nlog(n)), however, for convolution, it
1046 * is not the case. Convolving one image (a) to a kernel (b), can be done by
1047 * dividing image a to several blocks proportional to (b). Thus, one don't need
1048 * to do FFT for the whole image. The image can be divided to n/m part, and
1049 * the FFT itself is O(mlog(m)), so, the convolution process has time complexity
1050 * of O(nlog(m)) */
1051 if ((b->rows * b->cols < (log((double)(b->rows * b->cols)) + 1) * 15) && (a->type & CCV_8U))
1052 {
1053 _ccv_filter_direct_8u(a, b, dd, padding_pattern);
1054 } else {
1055#ifdef HAVE_FFTW31
1056 _ccv_filter_fftw(a, b, dd, padding_pattern);
1057#else
1058 _ccv_filter_kissfft(a, b, dd, padding_pattern);
1059#endif
1060 }
1061}
1062
1063void ccv_filter_kernel(ccv_dense_matrix_t* x, ccv_filter_kernel_f func, void* data)
1064{
1065 int i, j, k, ch = CCV_GET_CHANNEL(x->type)((x->type) & 0xFFF);
1066 unsigned char* m_ptr = x->data.u8;
1067 double rows_2 = (x->rows - 1) * 0.5;
1068 double cols_2 = (x->cols - 1) * 0.5;
1069#define for_block(_, _for_set) \
1070 for (i = 0; i < x->rows; i++) \
1071 { \
1072 for (j = 0; j < x->cols; j++) \
1073 { \
1074 double result = func(j - cols_2, i - rows_2, data); \
1075 for (k = 0; k < ch; k++) \
1076 _for_set(m_ptr, j * ch + k, result); \
1077 } \
1078 m_ptr += x->step; \
1079 }
1080 ccv_matrix_setter(x->type, for_block){ switch (((x->type) & 0xFF000)) { case CCV_32S: { for_block
(, _ccv_set_32s_value); break; } case CCV_32F: { for_block(, _ccv_set_32f_value
); break; } case CCV_64S: { for_block(, _ccv_set_64s_value); break
; } case CCV_64F: { for_block(, _ccv_set_64f_value); break; }
default: { for_block(, _ccv_set_8u_value); } } };
1081#undef for_block
1082 ccv_make_matrix_immutable(x);
1083}
1084
1085void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, ccv_dense_matrix_t** x, int x_type, ccv_dense_matrix_t** y, int y_type, double dx, double dy, double dxx, double dyy, int flag)
1086{
1087 assert(!(flag & CCV_L2_NORM) && (flag & CCV_GSEDT))((void) sizeof ((!(flag & CCV_L2_NORM) && (flag &
CCV_GSEDT)) ? 1 : 0), __extension__ ({ if (!(flag & CCV_L2_NORM
) && (flag & CCV_GSEDT)) ; else __assert_fail ("!(flag & CCV_L2_NORM) && (flag & CCV_GSEDT)"
, "ccv_numeric.c", 1087, __extension__ __PRETTY_FUNCTION__); }
));
1088 ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_distance_transform(%la,%la,%la,%la,%d)", dx, dy, dxx, dyy, flag), a->sig, CCV_EOF_SIGN)char _ccv_identifier_1088[(64)]; memset(_ccv_identifier_1088,
0, (64)); snprintf(_ccv_identifier_1088, (64), ("ccv_distance_transform(%la,%la,%la,%la,%d)"
), dx, dy, dxx, dyy, flag); size_t _ccv_string_size_1088 = (64
);; uint64_t sig = (a->sig != 0) ? ccv_cache_generate_signature
(_ccv_identifier_1088, _ccv_string_size_1088, a->sig, ((uint64_t
)0)) : 0;;
1089 type = (CCV_GET_DATA_TYPE(type)((type) & 0xFF000) == CCV_64F || CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) == CCV_64F || CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) == CCV_64S) ? CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF) | CCV_64F : CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF) | CCV_32F;
1090 ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_ALL_DATA_TYPE(CCV_8U | CCV_32S | CCV_32F | CCV_64S | CCV_64F | CCV_16F | CCV_QX
) | CCV_GET_CHANNEL(a->type)((a->type) & 0xFFF), type, sig);
1091 ccv_dense_matrix_t* mx = 0;
1092 ccv_dense_matrix_t* my = 0;
1093 if (x != 0)
1094 {
1095 ccv_declare_derived_signature(xsig, a->sig != 0, ccv_sign_with_format(64, "ccv_distance_transform_x(%la,%la,%la,%la,%d)", dx, dy, dxx, dyy, flag), a->sig, CCV_EOF_SIGN)char _ccv_identifier_1095[(64)]; memset(_ccv_identifier_1095,
0, (64)); snprintf(_ccv_identifier_1095, (64), ("ccv_distance_transform_x(%la,%la,%la,%la,%d)"
), dx, dy, dxx, dyy, flag); size_t _ccv_string_size_1095 = (64
);; uint64_t xsig = (a->sig != 0) ? ccv_cache_generate_signature
(_ccv_identifier_1095, _ccv_string_size_1095, a->sig, ((uint64_t
)0)) : 0;;
1096 mx = *x = ccv_dense_matrix_renew(*x, a->rows, a->cols, CCV_32S | CCV_C1, CCV_32S | CCV_C1, xsig);
1097 }
1098 if (y != 0)
1099 {
1100 ccv_declare_derived_signature(ysig, a->sig != 0, ccv_sign_with_format(64, "ccv_distance_transform_y(%la,%la,%la,%la,%d)", dx, dy, dxx, dyy, flag), a->sig, CCV_EOF_SIGN)char _ccv_identifier_1100[(64)]; memset(_ccv_identifier_1100,
0, (64)); snprintf(_ccv_identifier_1100, (64), ("ccv_distance_transform_y(%la,%la,%la,%la,%d)"
), dx, dy, dxx, dyy, flag); size_t _ccv_string_size_1100 = (64
);; uint64_t ysig = (a->sig != 0) ? ccv_cache_generate_signature
(_ccv_identifier_1100, _ccv_string_size_1100, a->sig, ((uint64_t
)0)) : 0;;
1101 my = *y = ccv_dense_matrix_renew(*y, a->rows, a->cols, CCV_32S | CCV_C1, CCV_32S | CCV_C1, ysig);
1102 }
1103 ccv_object_return_if_cached(, db, mx, my){ if ((!(db) || (((int*)(db))[0] & CCV_GARBAGE)) &&
(!(mx) || (((int*)(mx))[0] & CCV_GARBAGE)) && (!
(my) || (((int*)(my))[0] & CCV_GARBAGE)) && (!(0)
|| (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] &
CCV_GARBAGE)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE
)) && (!(0) || (((int*)(0))[0] & CCV_GARBAGE)) &&
(!(0) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0
) || (((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (
((int*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int
*)(0))[0] & CCV_GARBAGE)) && (!(0) || (((int*)(0)
)[0] & CCV_GARBAGE))) { (void)((db) && (((int*)(db
))[0] &= ~CCV_GARBAGE));(void)((mx) && (((int*)(mx
))[0] &= ~CCV_GARBAGE));(void)((my) && (((int*)(my
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));;; return ; } };
1104 ccv_revive_object_if_cached(db, mx, my)(void)((db) && (((int*)(db))[0] &= ~CCV_GARBAGE))
;(void)((mx) && (((int*)(mx))[0] &= ~CCV_GARBAGE)
);(void)((my) && (((int*)(my))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(
void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void
)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)
((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)(
(0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((
0) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0
) && (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0)
&& (((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) &&
(((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (
((int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && ((
(int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((
int*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int
*)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*
)(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)
(0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(
0))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0
))[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0)
)[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))
[0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[
0] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0
] &= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0]
&= ~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &=
~CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~
CCV_GARBAGE));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE
));(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE)
);(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE))
;(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));
(void)((0) && (((int*)(0))[0] &= ~CCV_GARBAGE));;;
1105 int i, j, k;
1106 unsigned char* a_ptr = a->data.u8;
1107 unsigned char* b_ptr = db->data.u8;
1108 int* v = (int*)alloca(sizeof(int) * ccv_max(db->rows, db->cols))__builtin_alloca (sizeof(int) * ({ typeof (db->rows) _a = (
db->rows); typeof (db->cols) _b = (db->cols); (_a >
_b) ? _a : _b; }));
1109 unsigned char* c_ptr = (unsigned char*)alloca(CCV_GET_DATA_TYPE_SIZE(db->type) * db->rows)__builtin_alloca (_ccv_get_data_type_size[((db->type) &
0xFF000) >> 12] * db->rows);
1110 int* x_ptr = mx ? mx->data.i32 : 0;
1111 int* y_ptr = my ? my->data.i32 : 0;
1112#define for_block(_for_max, _for_type_b, _for_set_b, _for_get_b, _for_get_a) \
1113 _for_type_b _dx = dx, _dy = dy, _dxx = dxx, _dyy = dyy; \
1114 _for_type_b* z = (_for_type_b*)alloca(sizeof(_for_type_b) * (ccv_max(db->rows, db->cols) + 1))__builtin_alloca (sizeof(_for_type_b) * (({ typeof (db->rows
) _a = (db->rows); typeof (db->cols) _b = (db->cols)
; (_a > _b) ? _a : _b; }) + 1)); \
1115 if (_dxx > 1e-6) \
1116 { \
1117 for (i = 0; i < a->rows; i++) \
1118 { \
1119 k = 0; \
1120 v[0] = 0; \
1121 z[0] = (_for_type_b)-_for_max; \
1122 z[1] = (_for_type_b)_for_max; \
1123 for (j = 1; j < a->cols; j++) \
1124 { \
1125 _for_type_b s; \
1126 for (;;) \
1127 { \
1128 assert(k >= 0 && k < ccv_max(db->rows, db->cols))((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; })) ? 1 : 0), __extension__ ({ if (
k >= 0 && k < ({ typeof (db->rows) _a = (db->
rows); typeof (db->cols) _b = (db->cols); (_a > _b) ?
_a : _b; })) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols)"
, "ccv_numeric.c", 1128, __extension__ __PRETTY_FUNCTION__); }
)); \
1129 s = ((SGN _for_get_a(a_ptr, j) + _dxx * j * j - _dx * j) - (SGN _for_get_a(a_ptr, v[k]) + _dxx * v[k] * v[k] - _dx * v[k])) / (2.0 * _dxx * (j - v[k])); \
1130 if (s > z[k]) break; \
1131 --k; \
1132 } \
1133 ++k; \
1134 assert(k >= 0 && k < ccv_max(db->rows, db->cols))((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; })) ? 1 : 0), __extension__ ({ if (
k >= 0 && k < ({ typeof (db->rows) _a = (db->
rows); typeof (db->cols) _b = (db->cols); (_a > _b) ?
_a : _b; })) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols)"
, "ccv_numeric.c", 1134, __extension__ __PRETTY_FUNCTION__); }
)); \
1135 v[k] = j; \
1136 z[k] = s; \
1137 z[k + 1] = (_for_type_b)_for_max; \
1138 } \
1139 assert(z[k + 1] >= a->cols - 1)((void) sizeof ((z[k + 1] >= a->cols - 1) ? 1 : 0), __extension__
({ if (z[k + 1] >= a->cols - 1) ; else __assert_fail (
"z[k + 1] >= a->cols - 1", "ccv_numeric.c", 1139, __extension__
__PRETTY_FUNCTION__); })); \
1140 k = 0; \
1141 if (mx) \
1142 { \
1143 for (j = 0; j < a->cols; j++) \
1144 { \
1145 while (z[k + 1] < j) \
1146 { \
1147 assert(k >= 0 && k < ccv_max(db->rows, db->cols) - 1)((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; }) - 1) ? 1 : 0), __extension__ ({
if (k >= 0 && k < ({ typeof (db->rows) _a =
(db->rows); typeof (db->cols) _b = (db->cols); (_a >
_b) ? _a : _b; }) - 1) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols) - 1"
, "ccv_numeric.c", 1147, __extension__ __PRETTY_FUNCTION__); }
)); \
1148 ++k; \
1149 } \
1150 _for_set_b(b_ptr, j, _dx * (j - v[k]) + _dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k])); \
1151 x_ptr[j] = j - v[k]; \
1152 } \
1153 x_ptr += mx->cols; \
1154 } else { \
1155 for (j = 0; j < a->cols; j++) \
1156 { \
1157 while (z[k + 1] < j) \
1158 { \
1159 assert(k >= 0 && k < ccv_max(db->rows, db->cols) - 1)((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; }) - 1) ? 1 : 0), __extension__ ({
if (k >= 0 && k < ({ typeof (db->rows) _a =
(db->rows); typeof (db->cols) _b = (db->cols); (_a >
_b) ? _a : _b; }) - 1) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols) - 1"
, "ccv_numeric.c", 1159, __extension__ __PRETTY_FUNCTION__); }
)); \
1160 ++k; \
1161 } \
1162 _for_set_b(b_ptr, j, _dx * (j - v[k]) + _dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k])); \
1163 } \
1164 } \
1165 a_ptr += a->step; \
1166 b_ptr += db->step; \
1167 } \
1168 } else { /* above algorithm cannot handle dxx == 0 properly, below is special casing for that */ \
1169 assert(mx == 0)((void) sizeof ((mx == 0) ? 1 : 0), __extension__ ({ if (mx ==
0) ; else __assert_fail ("mx == 0", "ccv_numeric.c", 1169, __extension__
__PRETTY_FUNCTION__); })); \
1170 for (i = 0; i < a->rows; i++) \
1171 { \
1172 for (j = 0; j < a->cols; j++) \
1173 _for_set_b(b_ptr, j, SGN _for_get_a(a_ptr, j)); \
1174 for (j = 1; j < a->cols; j++) \
1175 _for_set_b(b_ptr, j, ccv_min(_for_get_b(b_ptr, j), _for_get_b(b_ptr, j - 1) + _dx)({ typeof (_for_get_b(b_ptr, j)) _a = (_for_get_b(b_ptr, j));
typeof (_for_get_b(b_ptr, j - 1) + _dx) _b = (_for_get_b(b_ptr
, j - 1) + _dx); (_a < _b) ? _a : _b; })); \
1176 for (j = a->cols - 2; j >= 0; j--) \
1177 _for_set_b(b_ptr, j, ccv_min(_for_get_b(b_ptr, j), _for_get_b(b_ptr, j + 1) - _dx)({ typeof (_for_get_b(b_ptr, j)) _a = (_for_get_b(b_ptr, j));
typeof (_for_get_b(b_ptr, j + 1) - _dx) _b = (_for_get_b(b_ptr
, j + 1) - _dx); (_a < _b) ? _a : _b; })); \
1178 a_ptr += a->step; \
1179 b_ptr += db->step; \
1180 } \
1181 } \
1182 b_ptr = db->data.u8; \
1183 if (_dyy > 1e-6) \
1184 { \
1185 for (j = 0; j < db->cols; j++) \
1186 { \
1187 for (i = 0; i < db->rows; i++) \
1188 _for_set_b(c_ptr, i, _for_get_b(b_ptr + i * db->step, j)); \
1189 k = 0; \
1190 v[0] = 0; \
1191 z[0] = (_for_type_b)-_for_max; \
1192 z[1] = (_for_type_b)_for_max; \
1193 for (i = 1; i < db->rows; i++) \
1194 { \
1195 _for_type_b s; \
1196 for (;;) \
1197 { \
1198 assert(k >= 0 && k < ccv_max(db->rows, db->cols))((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; })) ? 1 : 0), __extension__ ({ if (
k >= 0 && k < ({ typeof (db->rows) _a = (db->
rows); typeof (db->cols) _b = (db->cols); (_a > _b) ?
_a : _b; })) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols)"
, "ccv_numeric.c", 1198, __extension__ __PRETTY_FUNCTION__); }
)); \
1199 s = ((_for_get_b(c_ptr, i) + _dyy * i * i - _dy * i) - (_for_get_b(c_ptr, v[k]) + _dyy * v[k] * v[k] - _dy * v[k])) / (2.0 * _dyy * (i - v[k])); \
1200 if (s > z[k]) break; \
1201 --k; \
1202 } \
1203 ++k; \
1204 assert(k >= 0 && k < ccv_max(db->rows, db->cols))((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; })) ? 1 : 0), __extension__ ({ if (
k >= 0 && k < ({ typeof (db->rows) _a = (db->
rows); typeof (db->cols) _b = (db->cols); (_a > _b) ?
_a : _b; })) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols)"
, "ccv_numeric.c", 1204, __extension__ __PRETTY_FUNCTION__); }
)); \
1205 v[k] = i; \
1206 z[k] = s; \
1207 z[k + 1] = (_for_type_b)_for_max; \
1208 } \
1209 assert(z[k + 1] >= db->rows - 1)((void) sizeof ((z[k + 1] >= db->rows - 1) ? 1 : 0), __extension__
({ if (z[k + 1] >= db->rows - 1) ; else __assert_fail (
"z[k + 1] >= db->rows - 1", "ccv_numeric.c", 1209, __extension__
__PRETTY_FUNCTION__); })); \
1210 k = 0; \
1211 if (my) \
1212 { \
1213 for (i = 0; i < db->rows; i++) \
1214 { \
1215 while (z[k + 1] < i) \
1216 { \
1217 assert(k >= 0 && k < ccv_max(db->rows, db->cols) - 1)((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; }) - 1) ? 1 : 0), __extension__ ({
if (k >= 0 && k < ({ typeof (db->rows) _a =
(db->rows); typeof (db->cols) _b = (db->cols); (_a >
_b) ? _a : _b; }) - 1) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols) - 1"
, "ccv_numeric.c", 1217, __extension__ __PRETTY_FUNCTION__); }
)); \
1218 ++k; \
1219 } \
1220 _for_set_b(b_ptr + i * db->step, j, _dy * (i - v[k]) + _dyy * (i - v[k]) * (i - v[k]) + _for_get_b(c_ptr, v[k])); \
1221 y_ptr[i * my->cols] = i - v[k]; \
1222 } \
1223 ++y_ptr; \
1224 } else { \
1225 for (i = 0; i < db->rows; i++) \
1226 { \
1227 while (z[k + 1] < i) \
1228 { \
1229 assert(k >= 0 && k < ccv_max(db->rows, db->cols) - 1)((void) sizeof ((k >= 0 && k < ({ typeof (db->
rows) _a = (db->rows); typeof (db->cols) _b = (db->cols
); (_a > _b) ? _a : _b; }) - 1) ? 1 : 0), __extension__ ({
if (k >= 0 && k < ({ typeof (db->rows) _a =
(db->rows); typeof (db->cols) _b = (db->cols); (_a >
_b) ? _a : _b; }) - 1) ; else __assert_fail ("k >= 0 && k < ccv_max(db->rows, db->cols) - 1"
, "ccv_numeric.c", 1229, __extension__ __PRETTY_FUNCTION__); }
)); \
1230 ++k; \
1231 } \
1232 _for_set_b(b_ptr + i * db->step, j, _dy * (i - v[k]) + _dyy * (i - v[k]) * (i - v[k]) + _for_get_b(c_ptr, v[k])); \
1233 } \
1234 } \
1235 } \
1236 } else { \
1237 assert(my == 0)((void) sizeof ((my == 0) ? 1 : 0), __extension__ ({ if (my ==
0) ; else __assert_fail ("my == 0", "ccv_numeric.c", 1237, __extension__
__PRETTY_FUNCTION__); })); \
1238 for (j = 0; j < db->cols; j++) \
1239 { \
1240 for (i = 1; i < db->rows; i++) \
1241 _for_set_b(b_ptr + i * db->step, j, ccv_min(_for_get_b(b_ptr + i * db->step, j), _for_get_b(b_ptr + (i - 1) * db->step, j) + _dy)({ typeof (_for_get_b(b_ptr + i * db->step, j)) _a = (_for_get_b
(b_ptr + i * db->step, j)); typeof (_for_get_b(b_ptr + (i -
1) * db->step, j) + _dy) _b = (_for_get_b(b_ptr + (i - 1)
* db->step, j) + _dy); (_a < _b) ? _a : _b; })); \
1242 for (i = db->rows - 2; i >= 0; i--) \
1243 _for_set_b(b_ptr + i * db->step, j, ccv_min(_for_get_b(b_ptr + i * db->step, j), _for_get_b(b_ptr + (i + 1) * db->step, j) - _dy)({ typeof (_for_get_b(b_ptr + i * db->step, j)) _a = (_for_get_b
(b_ptr + i * db->step, j)); typeof (_for_get_b(b_ptr + (i +
1) * db->step, j) - _dy) _b = (_for_get_b(b_ptr + (i + 1)
* db->step, j) - _dy); (_a < _b) ? _a : _b; })); \
1244 } \
1245 }
1246 if (flag & CCV_NEGATIVE)
1247 {
1248#define SGN -
1249 if (db->type & CCV_64F)
1250 {
1251 ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, DBL_MAX){ switch (((db->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, int
, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_8u_value);
} } }; break; } case CCV_32F: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, float
, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_8u_value);
} } }; break; } case CCV_64S: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_8u_value);
} } }; break; } case CCV_64F: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, double
, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_8u_value);
} } }; break; } default: { { switch (((a->type) & 0xFF000
)) { case CCV_32S: { for_block(1.7976931348623157e+308, unsigned
char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, unsigned
char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_8u_value
); } } }; } } };
1252 } else {
1253 ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, FLT_MAX){ switch (((db->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(3.40282347e+38F
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(3.40282347e+38F, int, _ccv_set_32s_value
, _ccv_get_32s_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(3.40282347e+38F, int, _ccv_set_32s_value, _ccv_get_32s_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(3.40282347e+38F
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64f_value
); break; } default: { for_block(3.40282347e+38F, int, _ccv_set_32s_value
, _ccv_get_32s_value, _ccv_get_8u_value); } } }; break; } case
CCV_32F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(3.40282347e+38F, float, _ccv_set_32f_value, _ccv_get_32f_value
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(3.40282347e+38F
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(3.40282347e+38F, float,
_ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64s_value);
break; } case CCV_64F: { for_block(3.40282347e+38F, float, _ccv_set_32f_value
, _ccv_get_32f_value, _ccv_get_64f_value); break; } default: {
for_block(3.40282347e+38F, float, _ccv_set_32f_value, _ccv_get_32f_value
, _ccv_get_8u_value); } } }; break; } case CCV_64S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(3.40282347e+38F
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32f_value)
; break; } case CCV_64S: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64s_value)
; break; } case CCV_64F: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64f_value)
; break; } default: { for_block(3.40282347e+38F, int64_t, _ccv_set_64s_value
, _ccv_get_64s_value, _ccv_get_8u_value); } } }; break; } case
CCV_64F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(3.40282347e+38F, double, _ccv_set_64f_value, _ccv_get_64f_value
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(3.40282347e+38F
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(3.40282347e+38F, double
, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64s_value)
; break; } case CCV_64F: { for_block(3.40282347e+38F, double,
_ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64f_value);
break; } default: { for_block(3.40282347e+38F, double, _ccv_set_64f_value
, _ccv_get_64f_value, _ccv_get_8u_value); } } }; break; } default
: { { switch (((a->type) & 0xFF000)) { case CCV_32S: {
for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value,
_ccv_get_8u_value, _ccv_get_32s_value); break; } case CCV_32F
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_64s_value); break; } case CCV_64F
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_64f_value); break; } default: {
for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value,
_ccv_get_8u_value, _ccv_get_8u_value); } } }; } } };
1254 }
1255#undef SGN
1256 } else {
1257#define SGN +
1258 if (db->type & CCV_64F)
1259 {
1260 ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, DBL_MAX){ switch (((db->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, int
, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_8u_value);
} } }; break; } case CCV_32F: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, float
, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_8u_value);
} } }; break; } case CCV_64S: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_8u_value);
} } }; break; } case CCV_64F: { { switch (((a->type) &
0xFF000)) { case CCV_32S: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, double
, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_8u_value);
} } }; break; } default: { { switch (((a->type) & 0xFF000
)) { case CCV_32S: { for_block(1.7976931348623157e+308, unsigned
char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_64s_value
); break; } case CCV_64F: { for_block(1.7976931348623157e+308
, unsigned char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_64f_value
); break; } default: { for_block(1.7976931348623157e+308, unsigned
char, _ccv_set_8u_value, _ccv_get_8u_value, _ccv_get_8u_value
); } } }; } } };
1261 } else {
1262 ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, FLT_MAX){ switch (((db->type) & 0xFF000)) { case CCV_32S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(3.40282347e+38F
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(3.40282347e+38F, int, _ccv_set_32s_value
, _ccv_get_32s_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(3.40282347e+38F, int, _ccv_set_32s_value, _ccv_get_32s_value
, _ccv_get_64s_value); break; } case CCV_64F: { for_block(3.40282347e+38F
, int, _ccv_set_32s_value, _ccv_get_32s_value, _ccv_get_64f_value
); break; } default: { for_block(3.40282347e+38F, int, _ccv_set_32s_value
, _ccv_get_32s_value, _ccv_get_8u_value); } } }; break; } case
CCV_32F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(3.40282347e+38F, float, _ccv_set_32f_value, _ccv_get_32f_value
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(3.40282347e+38F
, float, _ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(3.40282347e+38F, float,
_ccv_set_32f_value, _ccv_get_32f_value, _ccv_get_64s_value);
break; } case CCV_64F: { for_block(3.40282347e+38F, float, _ccv_set_32f_value
, _ccv_get_32f_value, _ccv_get_64f_value); break; } default: {
for_block(3.40282347e+38F, float, _ccv_set_32f_value, _ccv_get_32f_value
, _ccv_get_8u_value); } } }; break; } case CCV_64S: { { switch
(((a->type) & 0xFF000)) { case CCV_32S: { for_block(3.40282347e+38F
, int64_t, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32s_value
); break; } case CCV_32F: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_32f_value)
; break; } case CCV_64S: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64s_value)
; break; } case CCV_64F: { for_block(3.40282347e+38F, int64_t
, _ccv_set_64s_value, _ccv_get_64s_value, _ccv_get_64f_value)
; break; } default: { for_block(3.40282347e+38F, int64_t, _ccv_set_64s_value
, _ccv_get_64s_value, _ccv_get_8u_value); } } }; break; } case
CCV_64F: { { switch (((a->type) & 0xFF000)) { case CCV_32S
: { for_block(3.40282347e+38F, double, _ccv_set_64f_value, _ccv_get_64f_value
, _ccv_get_32s_value); break; } case CCV_32F: { for_block(3.40282347e+38F
, double, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_32f_value
); break; } case CCV_64S: { for_block(3.40282347e+38F, double
, _ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64s_value)
; break; } case CCV_64F: { for_block(3.40282347e+38F, double,
_ccv_set_64f_value, _ccv_get_64f_value, _ccv_get_64f_value);
break; } default: { for_block(3.40282347e+38F, double, _ccv_set_64f_value
, _ccv_get_64f_value, _ccv_get_8u_value); } } }; break; } default
: { { switch (((a->type) & 0xFF000)) { case CCV_32S: {
for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value,
_ccv_get_8u_value, _ccv_get_32s_value); break; } case CCV_32F
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_32f_value); break; } case CCV_64S
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_64s_value); break; } case CCV_64F
: { for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value
, _ccv_get_8u_value, _ccv_get_64f_value); break; } default: {
for_block(3.40282347e+38F, unsigned char, _ccv_set_8u_value,
_ccv_get_8u_value, _ccv_get_8u_value); } } }; } } };
1263 }
1264#undef SGN
1265 }
1266#undef for_block
1267}
1268
1269inline static double _kmeans1d_cost(double* cumsum, double* cumsum2, int i, int j)
1270{
1271 if (j < i)
1272 return 0;
1273 double mu = (cumsum[j + 1] - cumsum[i]) / (j - i + 1);
1274 double result = cumsum2[j + 1] - cumsum2[i];
1275 result += (j - i + 1) * (mu * mu);
1276 result -= (2 * mu) * (cumsum[j + 1] - cumsum[i]);
1277 return result;
1278}
1279
1280inline static double _kmeans1d_lookup(double* D, double* cumsum, double* cumsum2, int i, int j)
1281{
1282 const int col = i < j - 1 ? i : j - 1;
36
←
Assuming the condition is true
→
37
←
'?' condition is true
→
38
←
'col' initialized to 1
→
1283 return (col
38.1
'col' is >= 0
 >= 0 ? D[col] : 0) + _kmeans1d_cost(cumsum, cumsum2, j, i);
39
←
'?' condition is true
→
40
←
The left operand of '+' is a garbage value
1284}
1285
1286static void _smawk2(int row_start, int row_stride, int row_size, int* cols, int col_size, int* reserved, double* D, double* cumsum, double* cumsum2, int* result)
1287{
1288 if (row_size == 0)
25
←
Assuming 'row_size' is equal to 0
→
26
←
Taking true branch
→
1289 return;
1290 int* _cols = cols + col_size;
1291 int _col_size = 0;
1292 int i;
1293 for (i = 0; i < col_size; i++)
1294 {
1295 const int col = cols[i];
1296 for (;;)
1297 {
1298 if (_col_size == 0)
1299 break;
1300 const int row = row_start + row_stride * (_col_size - 1);
1301 if (_kmeans1d_lookup(D, cumsum, cumsum2, row, col) >= _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[_col_size - 1]))
1302 break;
1303 --_col_size;
1304 }
1305 if (_col_size < row_size)
1306 {
1307 _cols[_col_size] = col;
1308 ++_col_size;
1309 }
1310 }
1311 _smawk2(row_start + row_stride, row_stride * 2, row_size / 2, _cols, _col_size, reserved, D, cumsum, cumsum2, result);
1312 // Build the reverse lookup table.
1313 for (i = 0; i < _col_size; i++)
1314 reserved[_cols[i]] = i;
1315 int start = 0;
1316 for (i = 0; i < row_size; i += 2) {
1317 const int row = row_start + i * row_stride;
1318 int stop = _col_size - 1;
1319 if (i < row_size - 1)
1320 {
1321 const int argmin = result[row_start + (i + 1) * row_stride];
1322 stop = reserved[argmin];
1323 }
1324 int argmin = _cols[start];
1325 double min = _kmeans1d_lookup(D, cumsum, cumsum2, row, argmin);
1326 int c;
1327 for (c = start + 1; c <= stop; c++)
1328 {
1329 double value = _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[c]);
1330 if (value < min) {
1331 argmin = _cols[c];
1332 min = value;
1333 }
1334 }
1335 result[row] = argmin;
1336 start = stop;
1337 }
1338}
1339
1340static void _smawk1(int row_start, int row_stride, int row_size, int* cols, int col_size, int* reserved, double* D, double* cumsum, double* cumsum2, int* result)
1341{
1342 if (row_size == 0)
21
←
Assuming 'row_size' is not equal to 0
→
22
←
Taking false branch
→
1343 return;
1344 int* _cols = cols;
1345 int _col_size = 0;
1346 int i;
1347 for (i = 0; i
22.1
'i' is >= 'col_size'
 < col_size; i++)
23
←
Loop condition is false. Execution continues on line 1365
→
1348 {
1349 const int col = i;
1350 for (;;)
1351 {
1352 if (_col_size == 0)
1353 break;
1354 const int row = row_start + row_stride * (_col_size - 1);
1355 if (_kmeans1d_lookup(D, cumsum, cumsum2, row, col) >= _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[_col_size - 1]))
1356 break;
1357 --_col_size;
1358 }
1359 if (_col_size < row_size)
1360 {
1361 _cols[_col_size] = col;
1362 ++_col_size;
1363 }
1364 }
1365 _smawk2(row_start + row_stride, row_stride * 2, row_size / 2, _cols, _col_size, reserved, D, cumsum, cumsum2, result);
24
←
Calling '_smawk2'
→
27
←
Returning from '_smawk2'
→
1366 // Build the reverse lookup table.
1367 for (i = 0; i < _col_size; i++)
28
←
Loop condition is false. Execution continues on line 1369
→
1368 reserved[_cols[i]] = i;
1369 int start = 0;
1370 for (i = 0; i < row_size; i += 2) {
29
←
Assuming 'i' is < 'row_size'
→
30
←
Loop condition is true. Entering loop body
→
1371 const int row = row_start + i * row_stride;
31
←
'row' initialized to 1
→
1372 int stop = _col_size - 1;
1373 if (i < row_size - 1)
32
←
Assuming the condition is false
→
33
←
Taking false branch
→
1374 {
1375 const int argmin = result[row_start + (i + 1) * row_stride];
1376 stop = reserved[argmin];
1377 }
1378 int argmin = _cols[start];
1379 double min = _kmeans1d_lookup(D, cumsum, cumsum2, row, argmin);
34
←
Passing the value 1 via 4th parameter 'i'
→
35
←
Calling '_kmeans1d_lookup'
→
1380 int c;
1381 for (c = start + 1; c <= stop; c++)
1382 {
1383 double value = _kmeans1d_lookup(D, cumsum, cumsum2, row, _cols[c]);
1384 if (value < min) {
1385 argmin = _cols[c];
1386 min = value;
1387 }
1388 }
1389 result[row] = argmin;
1390 start = stop;
1391 }
1392}
1393
1394static void _smawk0(int row_start, int row_stride, int row_size, int* cols, int col_size, int* reserved, double* D, double* cumsum, double* cumsum2, int* result)
1395{
1396 if (row_size == 0)
18
←
Assuming 'row_size' is not equal to 0
→
19
←
Taking false branch
→
1397 return;
1398 _smawk1(row_start + row_stride, row_stride * 2, row_size / 2, cols, col_size, reserved, D, cumsum, cumsum2, result);
20
←
Calling '_smawk1'
→
1399 // Build the reverse lookup table.
1400 int start = 0;
1401 int i;
1402 for (i = 0; i < row_size; i += 2) {
1403 const int row = row_start + i * row_stride;
1404 int stop = col_size - 1;
1405 if (i < row_size - 1)
1406 {
1407 const int argmin = result[row_start + (i + 1) * row_stride];
1408 stop = argmin;
1409 }
1410 int argmin = start;
1411 double min = _kmeans1d_lookup(D, cumsum, cumsum2, row, argmin);
1412 int c;
1413 for (c = start + 1; c <= stop; c++)
1414 {
1415 double value = _kmeans1d_lookup(D, cumsum, cumsum2, row, c);
1416 if (value < min) {
1417 argmin = c;
1418 min = value;
1419 }
1420 }
1421 result[row] = argmin;
1422 start = stop;
1423 }
1424}
1425
1426static void smawk(const int n, int* cols, double* D, double* cumsum, double* cumsum2, int* result)
1427{
1428 _smawk0(0, 1, n, cols + n, n, cols, D, cumsum, cumsum2, result);
17
←
Calling '_smawk0'
→
1429}
1430
1431typedef struct {
1432 double value;
1433 int index;
1434} ccv_kmeans1d_undo_t;
1435
1436#undef more_than
1437#define less_than(s1, s2, aux) ((s1).value < (s2).value)
1438static CCV_IMPLEMENT_QSORT(_ccv_kmeans1d_undo_qsort, ccv_kmeans1d_undo_t, less_than)void _ccv_kmeans1d_undo_qsort(ccv_kmeans1d_undo_t *array, size_t
total, int aux) { int isort_thresh = 7; ccv_kmeans1d_undo_t t
; int sp = 0; struct { ccv_kmeans1d_undo_t *lb; ccv_kmeans1d_undo_t
*ub; } stack[48]; if( total <= 1 ) return; stack[0].lb = array
; stack[0].ub = array + (total - 1); while( sp >= 0 ) { ccv_kmeans1d_undo_t
* left = stack[sp].lb; ccv_kmeans1d_undo_t* right = stack[sp--
].ub; for(;;) { int i, n = (int)(right - left) + 1, m; ccv_kmeans1d_undo_t
* ptr; ccv_kmeans1d_undo_t* ptr2; if( n <= isort_thresh ) {
insert_sort: for( ptr = left + 1; ptr <= right; ptr++ ) {
for( ptr2 = ptr; ptr2 > left && less_than(ptr2[0]
,ptr2[-1], aux); ptr2--) (((t)) = ((ptr2[0])), ((ptr2[0])) = (
(ptr2[-1])), ((ptr2[-1])) = ((t))); } break; } else { ccv_kmeans1d_undo_t
* left0; ccv_kmeans1d_undo_t* left1; ccv_kmeans1d_undo_t* right0
; ccv_kmeans1d_undo_t* right1; ccv_kmeans1d_undo_t* pivot; ccv_kmeans1d_undo_t
* a; ccv_kmeans1d_undo_t* b; ccv_kmeans1d_undo_t* c; int swap_cnt
= 0; left0 = left; right0 = right; pivot = left + (n/2); if(
n > 40 ) { int d = n / 8; a = left, b = left + d, c = left
+ 2*d; left = less_than(*a, *b, aux) ? (less_than(*b, *c, aux
) ? b : (less_than(*a, *c, aux) ? c : a)) : (less_than(*c, *b
, aux) ? b : (less_than(*a, *c, aux) ? a : c)); a = pivot - d
, b = pivot, c = pivot + d; pivot = less_than(*a, *b, aux) ? (
less_than(*b, *c, aux) ? b : (less_than(*a, *c, aux) ? c : a)
) : (less_than(*c, *b, aux) ? b : (less_than(*a, *c, aux) ? a
: c)); a = right - 2*d, b = right - d, c = right; right = less_than
(*a, *b, aux) ? (less_than(*b, *c, aux) ? b : (less_than(*a, *
c, aux) ? c : a)) : (less_than(*c, *b, aux) ? b : (less_than(
*a, *c, aux) ? a : c)); } a = left, b = pivot, c = right; pivot
= less_than(*a, *b, aux) ? (less_than(*b, *c, aux) ? b : (less_than
(*a, *c, aux) ? c : a)) : (less_than(*c, *b, aux) ? b : (less_than
(*a, *c, aux) ? a : c)); if( pivot != left0 ) { (((t)) = ((*pivot
)), ((*pivot)) = ((*left0)), ((*left0)) = ((t))); pivot = left0
; } left = left1 = left0 + 1; right = right1 = right0; for(;;
) { while( left <= right && !less_than(*pivot, *left
, aux) ) { if( !less_than(*left, *pivot, aux) ) { if( left >
left1 ) (((t)) = ((*left1)), ((*left1)) = ((*left)), ((*left
)) = ((t))); swap_cnt = 1; left1++; } left++; } while( left <=
right && !less_than(*right, *pivot, aux) ) { if( !less_than
(*pivot, *right, aux) ) { if( right < right1 ) (((t)) = ((
*right1)), ((*right1)) = ((*right)), ((*right)) = ((t))); swap_cnt
= 1; right1--; } right--; } if( left > right ) break; (((
t)) = ((*left)), ((*left)) = ((*right)), ((*right)) = ((t)));
swap_cnt = 1; left++; right--; } if( swap_cnt == 0 ) { left =
left0, right = right0; goto insert_sort; } n = ({ typeof ((int
)(left1 - left0)) _a = ((int)(left1 - left0)); typeof ((int)(
left - left1)) _b = ((int)(left - left1)); (_a < _b) ? _a :
_b; }); for( i = 0; i < n; i++ ) (((t)) = ((left0[i])), (
(left0[i])) = ((left[i-n])), ((left[i-n])) = ((t))); n = ({ typeof
((int)(right0 - right1)) _a = ((int)(right0 - right1)); typeof
((int)(right1 - right)) _b = ((int)(right1 - right)); (_a <
_b) ? _a : _b; }); for( i = 0; i < n; i++ ) (((t)) = ((left
[i])), ((left[i])) = ((right0[i-n+1])), ((right0[i-n+1])) = (
(t))); n = (int)(left - left1); m = (int)(right1 - right); if
( n > 1 ) { if( m > 1 ) { if( n > m ) { stack[++sp].
lb = left0; stack[sp].ub = left0 + n - 1; left = right0 - m +
1, right = right0; } else { stack[++sp].lb = right0 - m + 1;
stack[sp].ub = right0; left = left0, right = left0 + n - 1; }
} else left = left0, right = left0 + n - 1; } else if( m >
1 ) left = right0 - m + 1, right = right0; else break; } } }
}
1439#undef less_than
1440
1441void ccv_kmeans1d(const ccv_dense_matrix_t* const a, const int k, int* const clusters, double* const centroids)
1442{
1443 assert(k > 1)((void) sizeof ((k > 1) ? 1 : 0), __extension__ ({ if (k >
1) ; else __assert_fail ("k > 1", "ccv_numeric.c", 1443, __extension__
__PRETTY_FUNCTION__); }));
1
Assuming 'k' is > 1
→
2
←
Taking true branch
→
1444 assert(CCV_GET_CHANNEL(a->type) == CCV_C1)((void) sizeof ((((a->type) & 0xFFF) == CCV_C1) ? 1 : 0
), __extension__ ({ if (((a->type) & 0xFFF) == CCV_C1)
; else __assert_fail ("CCV_GET_CHANNEL(a->type) == CCV_C1"
, "ccv_numeric.c", 1444, __extension__ __PRETTY_FUNCTION__); }
));
3
←
Assuming the condition is true
→
4
←
Taking true branch
→
1445 const int n = a->rows * a->cols;
1446 ccv_kmeans1d_undo_t* const sorted_undos = ccmallocmalloc(sizeof(ccv_kmeans1d_undo_t) * n);
1447 int i;
1448 if (CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) == CCV_16F)
5
←
Assuming the condition is false
→
6
←
Taking false branch
→
1449 {
1450 float* f = (float*)sorted_undos;
1451 ccv_half_precision_to_float((uint16_t*)a->data.f16, (float*)f, n);
1452 for (i = n - 1; i >= 0; i--)
1453 {
1454 sorted_undos[i].value = f[i];
1455 sorted_undos[i].index = i;
1456 }
1457 } else if (CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) == CCV_32F) {
7
←
Assuming the condition is false
→
8
←
Taking false branch
→
1458 for (i = 0; i < n; i++)
1459 {
1460 sorted_undos[i].value = a->data.f32[i];
1461 sorted_undos[i].index = i;
1462 }
1463 } else if (CCV_GET_DATA_TYPE(a->type)((a->type) & 0xFF000) == CCV_64F) {
9
←
Assuming the condition is false
→
10
←
Taking false branch
→
1464 for (i = 0; i < n; i++)
1465 {
1466 sorted_undos[i].value = a->data.f64[i];
1467 sorted_undos[i].index = i;
1468 }
1469 }
1470 _ccv_kmeans1d_undo_qsort(sorted_undos, n, 0);
1471 double* cc = ccmallocmalloc(sizeof(double) * 2 * (n + 1));
1472 cc[0] = 0;
1473 cc[n + 1] = 0;
1474 double cumsum = 0, cumsum2 = 0;
1475 for (i = 0; i < n; i++)
11
←
Assuming 'i' is >= 'n'
→
12
←
Loop condition is false. Execution continues on line 1482
→
1476 {
1477 cumsum += sorted_undos[i].value;
1478 cumsum2 += sorted_undos[i].value * sorted_undos[i].value;
1479 cc[i + 1] = cumsum;
1480 cc[i + 1 + n + 1] = cumsum2;
1481 }
1482 double* D = ccmallocmalloc(sizeof(double) * 2 * n);
13
←
Storing uninitialized value
→
1483 int* T = ccmallocmalloc(sizeof(int) * n * k);
1484 for (i = 0; i
13.1
'i' is >= 'n'
 < n; i++)
14
←
Loop condition is false. Execution continues on line 1489
→
1485 {
1486 D[i] = _kmeans1d_cost(cc, cc + n + 1, 0, i);
1487 T[i] = 0;
1488 }
1489 int k_;
1490 int* cols = ccmallocmalloc(sizeof(int) * n * 2);
1491 for (k_ = 1; k_
14.1
'k_' is < 'k'
 < k; k_++)
15
←
Loop condition is true. Entering loop body
→
1492 {
1493 double* lastD = D + ((k_ - 1) % 2) * n;
1494 int* const cT = T + k_ * n;
1495 smawk(n, cols, lastD, cc, cc + n + 1, cT);
16
←
Calling 'smawk'
→
1496 double* nextD = D + (k_ % 2) * n;
1497 for (i = 0; i < n; i++)
1498 {
1499 const int argmin = cT[i];
1500 nextD[i] = _kmeans1d_lookup(lastD, cc, cc + n + 1, i, argmin);
1501 }
1502 }
1503 ccfreefree(cc);
1504 ccfreefree(D);
1505 int t = n;
1506 k_ = k - 1;
1507 int n_ = n - 1;
1508 do {
1509 int t_ = t;
1510 t = T[k_ * n + n_];
1511 double centroid = 0.0;
1512 for (i = t; i < t_; i++)
1513 {
1514 cols[i] = k_;
1515 centroid += (sorted_undos[i].value - centroid) / (i - t + 1);
1516 }
1517 centroids[k_] = centroid;
1518 k_ -= 1;
1519 n_ = t - 1;
1520 } while (t > 0);
1521 for (i = 0; i < n; i++)
1522 clusters[sorted_undos[i].index] = cols[i];
1523 ccfreefree(cols);
1524 ccfreefree(sorted_undos);
1525 ccfreefree(T);
1526}