| File: | ccv_numeric.c |
| Warning: | line 1251, column 4 The left operand of '<' is a garbage value |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 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 | ||||
| 12 | const 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 | ||||
| 21 | void ccv_invert(ccv_matrix_t* a, ccv_matrix_t** b, int type) | |||
| 22 | { | |||
| 23 | } | |||
| 24 | ||||
| 25 | void ccv_solve(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type) | |||
| 26 | { | |||
| 27 | } | |||
| 28 | ||||
| 29 | void 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 | ||||
| 130 | void 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 */ | |||
| 386 | static 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 | ||||
| 567 | static 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 | ||||
| 583 | static pthread_mutex_t fftw_plan_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, PTHREAD_MUTEX_TIMED_NP, 0, 0, { 0, 0 } } }; | |||
| 584 | ||||
| 585 | static 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 | |||
| 771 | static 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 | ||||
| 960 | static 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 | ||||
| 1036 | void 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 | ||||
| 1063 | void 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 | ||||
| 1085 | void 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
| |||
| 1111 | int* y_ptr = my
| |||
| 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 | ||||
| 1269 | inline 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 | ||||
| 1280 | inline 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; | |||
| 1283 | return (col >= 0 ? D[col] : 0) + _kmeans1d_cost(cumsum, cumsum2, j, i); | |||
| 1284 | } | |||
| 1285 | ||||
| 1286 | static 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) | |||
| 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 | ||||
| 1340 | static 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) | |||
| 1343 | return; | |||
| 1344 | int* _cols = cols; | |||
| 1345 | int _col_size = 0; | |||
| 1346 | int i; | |||
| 1347 | for (i = 0; i < col_size; i++) | |||
| 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); | |||
| 1366 | // Build the reverse lookup table. | |||
| 1367 | for (i = 0; i < _col_size; i++) | |||
| 1368 | reserved[_cols[i]] = i; | |||
| 1369 | int start = 0; | |||
| 1370 | for (i = 0; i < row_size; i += 2) { | |||
| 1371 | const int row = row_start + i * row_stride; | |||
| 1372 | int stop = _col_size - 1; | |||
| 1373 | if (i < row_size - 1) | |||
| 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); | |||
| 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 | ||||
| 1394 | static 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) | |||
| 1397 | return; | |||
| 1398 | _smawk1(row_start + row_stride, row_stride * 2, row_size / 2, cols, col_size, reserved, D, cumsum, cumsum2, result); | |||
| 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 | ||||
| 1426 | static 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); | |||
| 1429 | } | |||
| 1430 | ||||
| 1431 | typedef 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) | |||
| 1438 | static 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 | ||||
| 1441 | void 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__); })); | |||
| 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__); } )); | |||
| 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) | |||
| 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) { | |||
| 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) { | |||
| 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++) | |||
| 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); | |||
| 1483 | int* T = ccmallocmalloc(sizeof(int) * n * k); | |||
| 1484 | for (i = 0; i < n; i++) | |||
| 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_ < k; k_++) | |||
| 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); | |||
| 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 | } |