#include "ccv.h"
#include "ccv_internal.h"
#if defined(HAVE_SSE2)
#include <xmmintrin.h>
#elif defined(HAVE_NEON)
#include <arm_neon.h>
#endif
/* sobel filter is fundamental to many other high-level algorithms,
* here includes 2 special case impl (for 1x3/3x1, 3x3) and one general impl */
void ccv_sobel(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int dx, int dy)
{
ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_sobel(%d,%d)", dx, dy), a->sig, CCV_EOF_SIGN);
type = (type == 0) ? CCV_32S | CCV_GET_CHANNEL(a->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(a->type);
ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_GET_CHANNEL(a->type) | CCV_ALL_DATA_TYPE, type, sig);
ccv_object_return_if_cached(, db);
int i, j, k, c, ch = CCV_GET_CHANNEL(a->type);
unsigned char* a_ptr = a->data.u8;
unsigned char* b_ptr = db->data.u8;
if (dx == 1 && dy == 0)
{
assert(a->cols >= 3);
/* special case 1: 1x3 or 3x1 window */
#define for_block(_for_get, _for_set) \
for (i = 0; i < a->rows; i++) \
{ \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, k, 2 * (_for_get(a_ptr, ch + k, 0) - _for_get(a_ptr, k, 0)), 0); \
for (j = 1; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, _for_get(a_ptr, (j + 1) * ch + k, 0) - _for_get(a_ptr, (j - 1) * ch + k, 0), 0); \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k, 0) - _for_get(a_ptr, (a->cols - 2) * ch + k, 0)), 0); \
b_ptr += db->step; \
a_ptr += a->step; \
}
ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);
#undef for_block
} else if (dx == 0 && dy == 1) {
assert(a->rows >= 3);
/* special case 1: 1x3 or 3x1 window */
#define for_block(_for_get, _for_set) \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, j * ch + k, 0) - _for_get(a_ptr, j * ch + k, 0)), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
for (i = 1; i < a->rows - 1; i++) \
{ \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k, 0) - _for_get(a_ptr - a->step, j * ch + k, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
} \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k, 0) - _for_get(a_ptr - a->step, j * ch + k, 0)), 0);
ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);
#undef for_block
} else if ((dx == 1 && dy == 1) || (dx == -1 && dy == -1)) {
/* special case 2: 3x3 window with diagonal direction */
assert(a->rows >= 3 && a->cols >= 3);
#define for_block(_for_get, _for_set) \
for (j = 0; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, (j + 1) * ch + k, 0) - _for_get(a_ptr, j * ch + k, 0)), 0); \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr + a->step, (a->cols - 1) * ch + k, 0) - _for_get(a_ptr, (a->cols - 1) * ch + k, 0)), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
for (i = 1; i < a->rows - 1; i++) \
{ \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, k, 2 * (_for_get(a_ptr + a->step, ch + k, 0) - _for_get(a_ptr, k, 0)), 0); \
for (j = 1; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, _for_get(a_ptr + a->step, (j + 1) * ch + k, 0) - _for_get(a_ptr - a->step, (j - 1) * ch + k, 0), 0); \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k, 0) - _for_get(a_ptr - a->step, (a->cols - 2) * ch + k, 0)), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
} \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, k, 2 * (_for_get(a_ptr, k, 0) - _for_get(a_ptr - a->step, k, 0)), 0); \
for (j = 1; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k, 0) - _for_get(a_ptr - a->step, (j - 1) * ch + k, 0)), 0);
ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);
#undef for_block
} else if ((dx == 1 && dy == -1) || (dx == -1 && dy == 1)) {
/* special case 2: 3x3 window with diagonal direction */
assert(a->rows >= 3 && a->cols >= 3);
#define for_block(_for_get, _for_set) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, k, 2 * (_for_get(a_ptr + a->step, k, 0) - _for_get(a_ptr, k, 0)), 0); \
for (j = 1; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, (j - 1) * ch + k, 0) - _for_get(a_ptr, j * ch + k, 0)), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
for (i = 1; i < a->rows - 1; i++) \
{ \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, k, 2 * (_for_get(a_ptr, k, 0) - _for_get(a_ptr - a->step, ch + k, 0)), 0); \
for (j = 1; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, _for_get(a_ptr + a->step, (j - 1) * ch + k, 0) - _for_get(a_ptr - a->step, (j + 1) * ch + k, 0), 0); \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr + a->step, (a->cols - 2) * ch + k, 0) - _for_get(a_ptr, (a->cols - 1) * ch + k, 0)), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
} \
for (j = 0; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k, 0) - _for_get(a_ptr - a->step, (j + 1) * ch + k, 0)), 0); \
for (k = 0; k < ch; k++) \
_for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k, 0) - _for_get(a_ptr - a->step, (a->cols - 1) * ch + k, 0)), 0);
ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);
#undef for_block
} else if (dx == 3 && dy == 0) {
assert(a->rows >= 3 && a->cols >= 3);
/* special case 3: 3x3 window, corresponding sigma = 0.85 */
unsigned char* buf = (unsigned char*)alloca(db->step);
#define for_block(_for_get, _for_set_b, _for_get_b) \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k, 0) + 3 * _for_get(a_ptr, j * ch + k, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
for (i = 1; i < a->rows - 1; i++) \
{ \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k, 0) + 2 * _for_get(a_ptr, j * ch + k, 0) + _for_get(a_ptr - a->step, j * ch + k, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
} \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, 3 * _for_get(a_ptr, j * ch + k, 0) + _for_get(a_ptr - a->step, j * ch + k, 0), 0); \
b_ptr = db->data.u8; \
for (i = 0; i < a->rows; i++) \
{ \
for (k = 0; k < ch; k++) \
_for_set_b(buf, k, _for_get_b(b_ptr, ch + k, 0) - _for_get_b(b_ptr, k, 0), 0); \
for (j = 1; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + k, _for_get_b(b_ptr, (j + 1) * ch + k, 0) - _for_get_b(b_ptr, (j - 1) * ch + k, 0), 0); \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (a->cols - 1) * ch + k, _for_get_b(b_ptr, (a->cols - 1) * ch + k, 0) - _for_get_b(b_ptr, (a->cols - 2) * ch + k, 0), 0); \
memcpy(b_ptr, buf, db->step); \
b_ptr += db->step; \
}
ccv_matrix_getter(a->type, ccv_matrix_setter_getter, db->type, for_block);
#undef for_block
} else if (dx == 0 && dy == 3) {
assert(a->rows >= 3 && a->cols >= 3);
/* special case 3: 3x3 window, corresponding sigma = 0.85 */
unsigned char* buf = (unsigned char*)alloca(db->step);
#define for_block(_for_get, _for_set_b, _for_get_b) \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k, 0) - _for_get(a_ptr, j * ch + k, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
for (i = 1; i < a->rows - 1; i++) \
{ \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k, 0) - _for_get(a_ptr - a->step, j * ch + k, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
} \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(b_ptr, j * ch + k, _for_get(a_ptr, j * ch + k, 0) - _for_get(a_ptr - a->step, j * ch + k, 0), 0); \
b_ptr = db->data.u8; \
for (i = 0; i < a->rows; i++) \
{ \
for (k = 0; k < ch; k++) \
_for_set_b(buf, k, _for_get_b(b_ptr, ch + k, 0) + 3 * _for_get_b(b_ptr, k, 0), 0); \
for (j = 1; j < a->cols - 1; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + k, _for_get_b(b_ptr, (j + 1) * ch + k, 0) + 2 * _for_get_b(b_ptr, j * ch + k, 0) + _for_get_b(b_ptr, (j - 1) * ch + k, 0), 0); \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (a->cols - 1) * ch + k, _for_get_b(b_ptr, (a->cols - 2) * ch + k, 0) + 3 * _for_get_b(b_ptr, (a->cols - 1) * ch + k, 0), 0); \
memcpy(b_ptr, buf, db->step); \
b_ptr += db->step; \
}
ccv_matrix_getter(a->type, ccv_matrix_setter_getter, db->type, for_block);
#undef for_block
} else {
/* general case: in this case, I will generate a separable filter, and do the convolution */
int fsz = ccv_max(dx, dy);
assert(fsz % 2 == 1);
int hfz = fsz / 2;
unsigned char* df = (unsigned char*)alloca(sizeof(double) * fsz);
unsigned char* gf = (unsigned char*)alloca(sizeof(double) * fsz);
/* the sigma calculation is linear derviation of 3x3 - 0.85, 5x5 - 1.32 */
double sigma = ((fsz - 1) / 2) * 0.47 + 0.38;
double sigma2 = (2.0 * sigma * sigma);
/* 2.5 is the factor to make the kernel "visible" in integer setting */
double psigma3 = 2.5 / sqrt(sqrt(2 * CCV_PI) * sigma * sigma * sigma);
for (i = 0; i < fsz; i++)
{
((double*)df)[i] = (i - hfz) * exp(-((i - hfz) * (i - hfz)) / sigma2) * psigma3;
((double*)gf)[i] = exp(-((i - hfz) * (i - hfz)) / sigma2) * psigma3;
}
if (db->type & CCV_32S)
{
for (i = 0; i < fsz; i++)
{
// df could be negative, thus, (int)(x + 0.5) shortcut will not work
((int*)df)[i] = (int)round(((double*)df)[i] * 256.0);
((int*)gf)[i] = (int)(((double*)gf)[i] * 256.0 + 0.5);
}
} else {
for (i = 0; i < fsz; i++)
{
ccv_set_value(db->type, df, i, ((double*)df)[i], 0);
ccv_set_value(db->type, gf, i, ((double*)gf)[i], 0);
}
}
if (dx < dy)
{
unsigned char* tf = df;
df = gf;
gf = tf;
}
unsigned char* buf = (unsigned char*)alloca(sizeof(double) * ch * (fsz + ccv_max(a->rows, a->cols)));
#define for_block(_for_get, _for_type_b, _for_set_b, _for_get_b) \
for (i = 0; i < a->rows; i++) \
{ \
for (j = 0; j < hfz; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + k, _for_get(a_ptr, k, 0), 0); \
for (j = 0; j < a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (j + hfz) * ch + k, _for_get(a_ptr, j * ch + k, 0), 0); \
for (j = a->cols; j < a->cols + hfz; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (j + hfz) * ch + k, _for_get(a_ptr, (a->cols - 1) * ch + k, 0), 0); \
for (j = 0; j < a->cols; j++) \
{ \
for (c = 0; c < ch; c++) \
{ \
_for_type_b sum = 0; \
for (k = 0; k < fsz; k++) \
sum += _for_get_b(buf, (j + k) * ch + c, 0) * _for_get_b(df, k, 0); \
_for_set_b(b_ptr, j * ch + c, sum, 8); \
} \
} \
a_ptr += a->step; \
b_ptr += db->step; \
} \
b_ptr = db->data.u8; \
for (i = 0; i < a->cols; i++) \
{ \
for (j = 0; j < hfz; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + k, _for_get_b(b_ptr, i * ch + k, 0), 0); \
for (j = 0; j < a->rows; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + j * db->step, i * ch + k, 0), 0); \
for (j = a->rows; j < a->rows + hfz; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + (a->rows - 1) * db->step, i * ch + k, 0), 0); \
for (j = 0; j < a->rows; j++) \
{ \
for (c = 0; c < ch; c++) \
{ \
_for_type_b sum = 0; \
for (k = 0; k < fsz; k++) \
sum += _for_get_b(buf, (j + k) * ch + c, 0) * _for_get_b(gf, k, 0); \
_for_set_b(b_ptr + j * db->step, i * ch + c, sum, 8); \
} \
} \
}
ccv_matrix_getter(a->type, ccv_matrix_typeof_setter_getter, db->type, for_block);
#undef for_block
}
}
/* the fast arctan function adopted from OpenCV */
static void _ccv_atan2(float* x, float* y, float* angle, float* mag, int len)
{
int i = 0;
float scale = (float)(180.0 / CCV_PI);
#ifdef HAVE_SSE2
#ifndef _WIN32
union { int i; float fl; } iabsmask; iabsmask.i = 0x7fffffff;
__m128 eps = _mm_set1_ps((float)1e-6), absmask = _mm_set1_ps(iabsmask.fl);
__m128 _90 = _mm_set1_ps((float)(3.141592654 * 0.5)), _180 = _mm_set1_ps((float)3.141592654), _360 = _mm_set1_ps((float)(3.141592654 * 2));
__m128 zero = _mm_setzero_ps(), _0_28 = _mm_set1_ps(0.28f), scale4 = _mm_set1_ps(scale);
for(; i <= len - 4; i += 4)
{
__m128 x4 = _mm_loadu_ps(x + i), y4 = _mm_loadu_ps(y + i);
__m128 xq4 = _mm_mul_ps(x4, x4), yq4 = _mm_mul_ps(y4, y4);
__m128 xly = _mm_cmplt_ps(xq4, yq4);
__m128 z4 = _mm_div_ps(_mm_mul_ps(x4, y4), _mm_add_ps(_mm_add_ps(_mm_max_ps(xq4, yq4), _mm_mul_ps(_mm_min_ps(xq4, yq4), _0_28)), eps));
// a4 <- x < y ? 90 : 0;
__m128 a4 = _mm_and_ps(xly, _90);
// a4 <- (y < 0 ? 360 - a4 : a4) == ((x < y ? y < 0 ? 270 : 90) : (y < 0 ? 360 : 0))
__m128 mask = _mm_cmplt_ps(y4, zero);
a4 = _mm_or_ps(_mm_and_ps(_mm_sub_ps(_360, a4), mask), _mm_andnot_ps(mask, a4));
// a4 <- (x < 0 && !(x < y) ? 180 : a4)
mask = _mm_andnot_ps(xly, _mm_cmplt_ps(x4, zero));
a4 = _mm_or_ps(_mm_and_ps(_180, mask), _mm_andnot_ps(mask, a4));
// a4 <- (x < y ? a4 - z4 : a4 + z4)
a4 = _mm_mul_ps(_mm_add_ps(_mm_xor_ps(z4, _mm_andnot_ps(absmask, xly)), a4), scale4);
__m128 m4 = _mm_sqrt_ps(_mm_add_ps(xq4, yq4));
_mm_storeu_ps(angle + i, a4);
_mm_storeu_ps(mag + i, m4);
}
#endif
#endif
for(; i < len; i++)
{
float xf = x[i], yf = y[i];
float a, x2 = xf * xf, y2 = yf * yf;
if(y2 <= x2)
a = xf * yf / (x2 + 0.28f * y2 + (float)1e-6) + (float)(xf < 0 ? CCV_PI : yf >= 0 ? 0 : CCV_PI * 2);
else
a = (float)(yf >= 0 ? CCV_PI * 0.5 : CCV_PI * 1.5) - xf * yf / (y2 + 0.28f * x2 + (float)1e-6);
angle[i] = a * scale;
mag[i] = sqrtf(x2 + y2);
}
}
void ccv_gradient(ccv_dense_matrix_t* a, ccv_dense_matrix_t** theta, int ttype, ccv_dense_matrix_t** m, int mtype, int dx, int dy)
{
ccv_declare_derived_signature(tsig, a->sig != 0, ccv_sign_with_format(64, "ccv_gradient(theta,%d,%d)", dx, dy), a->sig, CCV_EOF_SIGN);
ccv_declare_derived_signature(msig, a->sig != 0, ccv_sign_with_format(64, "ccv_gradient(m,%d,%d)", dx, dy), a->sig, CCV_EOF_SIGN);
int ch = CCV_GET_CHANNEL(a->type);
ccv_dense_matrix_t* dtheta = *theta = ccv_dense_matrix_renew(*theta, a->rows, a->cols, CCV_32F | ch, CCV_32F | ch, tsig);
ccv_dense_matrix_t* dm = *m = ccv_dense_matrix_renew(*m, a->rows, a->cols, CCV_32F | ch, CCV_32F | ch, msig);
assert(dtheta && dm);
ccv_object_return_if_cached(, dtheta, dm);
ccv_revive_object_if_cached(dtheta, dm);
ccv_dense_matrix_t* tx = 0;
ccv_dense_matrix_t* ty = 0;
ccv_sobel(a, &tx, CCV_32F | ch, dx, 0);
ccv_sobel(a, &ty, CCV_32F | ch, 0, dy);
_ccv_atan2(tx->data.f32, ty->data.f32, dtheta->data.f32, dm->data.f32, ch * a->rows * a->cols);
ccv_matrix_free(tx);
ccv_matrix_free(ty);
}
static void _ccv_flip_y_self(ccv_dense_matrix_t* a)
{
int i;
unsigned char* buffer = (unsigned char*)alloca(a->step);
unsigned char* a_ptr = a->data.u8;
unsigned char* b_ptr = a->data.u8 + (a->rows - 1) * a->step;
for (i = 0; i < a->rows / 2; i++)
{
memcpy(buffer, a_ptr, a->step);
memcpy(a_ptr, b_ptr, a->step);
memcpy(b_ptr, buffer, a->step);
a_ptr += a->step;
b_ptr -= a->step;
}
}
static void _ccv_flip_x_self(ccv_dense_matrix_t* a)
{
int i, j;
int len = CCV_GET_DATA_TYPE_SIZE(a->type) * CCV_GET_CHANNEL(a->type);
unsigned char* buffer = (unsigned char*)alloca(len);
unsigned char* a_ptr = a->data.u8;
for (i = 0; i < a->rows; i++)
{
for (j = 0; j < a->cols / 2; j++)
{
memcpy(buffer, a_ptr + j * len, len);
memcpy(a_ptr + j * len, a_ptr + (a->cols - 1 - j) * len, len);
memcpy(a_ptr + (a->cols - 1 - j) * len, buffer, len);
}
a_ptr += a->step;
}
}
void ccv_flip(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int btype, int type)
{
/* this is the special case where ccv_declare_derived_signature_* macros cannot handle properly */
uint64_t sig = a->sig;
if (type & CCV_FLIP_Y)
sig = (a->sig == 0) ? 0 : ccv_cache_generate_signature("ccv_flip_y", 10, sig, CCV_EOF_SIGN);
if (type & CCV_FLIP_X)
sig = (a->sig == 0) ? 0 : ccv_cache_generate_signature("ccv_flip_x", 10, sig, CCV_EOF_SIGN);
ccv_dense_matrix_t* db;
if (b == 0)
{
db = a;
if (a->sig != 0)
{
btype = CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type);
sig = ccv_cache_generate_signature((const char*)&btype, sizeof(int), sig, CCV_EOF_SIGN);
a->sig = sig;
}
} else {
btype = CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type);
*b = db = ccv_dense_matrix_renew(*b, a->rows, a->cols, btype, btype, sig);
ccv_object_return_if_cached(, db);
if (a->data.u8 != db->data.u8)
memcpy(db->data.u8, a->data.u8, a->rows * a->step);
}
if (type & CCV_FLIP_Y)
_ccv_flip_y_self(db);
if (type & CCV_FLIP_X)
_ccv_flip_x_self(db);
}
void ccv_blur(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, double sigma)
{
ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_blur(%la)", sigma), a->sig, CCV_EOF_SIGN);
type = (type == 0) ? CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(a->type);
ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(a->type), type, sig);
ccv_object_return_if_cached(, db);
int fsz = ccv_max(1, (int)(4.0 * sigma + 1.0 - 1e-8)) * 2 + 1;
int hfz = fsz / 2;
assert(hfz > 0);
unsigned char* buf = (unsigned char*)alloca(sizeof(double) * ccv_max(hfz * 2 + a->rows, (hfz * 2 + a->cols) * CCV_GET_CHANNEL(a->type)));
unsigned char* filter = (unsigned char*)alloca(sizeof(double) * fsz);
double tw = 0;
int i, j, k, ch = CCV_GET_CHANNEL(a->type);
assert(fsz > 0);
for (i = 0; i < fsz; i++)
tw += ((double*)filter)[i] = exp(-((i - hfz) * (i - hfz)) / (2.0 * sigma * sigma));
int no_8u_type = (db->type & CCV_8U) ? CCV_32S : db->type;
if (no_8u_type & CCV_32S)
{
tw = 256.0 / tw;
for (i = 0; i < fsz; i++)
((int*)filter)[i] = (int)(((double*)filter)[i] * tw + 0.5);
} else {
tw = 1.0 / tw;
for (i = 0; i < fsz; i++)
ccv_set_value(no_8u_type, filter, i, ((double*)filter)[i] * tw, 0);
}
/* horizontal */
unsigned char* a_ptr = a->data.u8;
unsigned char* b_ptr = db->data.u8;
assert(ch > 0);
#define for_block(_for_type, _for_set_b, _for_get_b, _for_set_a, _for_get_a) \
for (i = 0; i < a->rows; i++) \
{ \
for (j = 0; j < hfz; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + k, _for_get_a(a_ptr, k, 0), 0); \
for (j = 0; j < a->cols * ch; j++) \
_for_set_b(buf, j + hfz * ch, _for_get_a(a_ptr, j, 0), 0); \
for (j = a->cols; j < hfz + a->cols; j++) \
for (k = 0; k < ch; k++) \
_for_set_b(buf, j * ch + hfz * ch + k, _for_get_a(a_ptr, (a->cols - 1) * ch + k, 0), 0); \
for (j = 0; j < a->cols * ch; j++) \
{ \
_for_type sum = 0; \
for (k = 0; k < fsz; k++) \
sum += _for_get_b(buf, k * ch + j, 0) * _for_get_b(filter, k, 0); \
_for_set_b(buf, j, sum, 8); \
} \
for (j = 0; j < a->cols * ch; j++) \
_for_set_a(b_ptr, j, _for_get_b(buf, j, 0), 0); \
a_ptr += a->step; \
b_ptr += db->step; \
}
ccv_matrix_typeof_setter_getter(no_8u_type, ccv_matrix_setter, db->type, ccv_matrix_getter, a->type, for_block);
#undef for_block
/* vertical */
b_ptr = db->data.u8;
#define for_block(_for_type, _for_set_b, _for_get_b, _for_set_a, _for_get_a) \
for (i = 0; i < a->cols * ch; i++) \
{ \
for (j = 0; j < hfz; j++) \
_for_set_b(buf, j, _for_get_a(b_ptr, i, 0), 0); \
for (j = 0; j < a->rows; j++) \
_for_set_b(buf, j + hfz, _for_get_a(b_ptr + j * db->step, i, 0), 0); \
for (j = a->rows; j < hfz + a->rows; j++) \
_for_set_b(buf, j + hfz, _for_get_a(b_ptr + (a->rows - 1) * db->step, i, 0), 0); \
for (j = 0; j < a->rows; j++) \
{ \
_for_type sum = 0; \
for (k = 0; k < fsz; k++) \
sum += _for_get_b(buf, k + j, 0) * _for_get_b(filter, k, 0); \
_for_set_b(buf, j, sum, 8); \
} \
for (j = 0; j < a->rows; j++) \
_for_set_a(b_ptr + j * db->step, i, _for_get_b(buf, j, 0), 0); \
}
ccv_matrix_typeof_setter_getter(no_8u_type, ccv_matrix_setter_getter, db->type, for_block);
#undef for_block
}