#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_TYPE0
(type) | 0
CCV_GET_CHANNEL0
(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 == 02
)

  {

    assert(a->cols >= 3);

    /* special case 1: 1x3 or 3x1 window */

#define for_block(_for_get, _for_set) \

    
for (i = 0; 1
i < a->rows; 
i++500
) \

    { \

      for (k = 0; k < ch; 
k++500
) \

        _for_set(b_ptr, k, 2 * (_for_get(a_ptr, ch + k) - _for_get(a_ptr, k))); \

      for (j = 1; j < a->cols - 1; 
j++249k
) \

        
for (k = 0; 249k
k < ch; 
k++249k
) \

          
_for_set500
(b_ptr, j * ch + k, _for_get(a_ptr, (j + 1) * ch + k) - _for_get(a_ptr, (j - 1) * ch + k)); \

      for (k = 0; k < ch; 
k++500
) \

        _for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k) - _for_get(a_ptr, (a->cols - 2) * ch + k))); \

      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 == 14
) {

    assert(a->rows >= 3);

    /* special case 1: 1x3 or 3x1 window */

#define for_block(_for_get, _for_set) \

    
for (j = 0; 1
j < a->cols; 
j++500
) \

      
for (k = 0; 500
k < ch; 
k++500
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, j * ch + k) - _for_get(a_ptr, j * ch + k))); \

    a_ptr += a->step; \

    b_ptr += db->step; \

    for (i = 1; i < a->rows - 1; 
i++498
) \

    { \

      for (j = 0; j < a->cols; 
j++249k
) \

        
for (k = 0; 249k
k < ch; 
k++249k
) \

          
_for_set498
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k) - _for_get(a_ptr - a->step, j * ch + k)); \

      a_ptr += a->step; \

      b_ptr += db->step; \

    } \

    for (j = 0; j < a->cols; 
j++500
) \

      
for (k = 0; 500
k < ch; 
k++500
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k) - _for_get(a_ptr - a->step, j * ch + k)));

    ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);

#undef for_block

  } else if ((dx == 1 && 
dy == 11
) || 
(7
dx == -17
 && 
dy == -11
)) {

    /* 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; 1
j < a->cols - 1; 
j++499
) \

      
for (k = 0; 499
k < ch; 
k++499
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, (j + 1) * ch + k) - _for_get(a_ptr, j * ch + k))); \

    for (k = 0; k < ch; 
k++1
) \

      _for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr + a->step, (a->cols - 1) * ch + k) - _for_get(a_ptr, (a->cols - 1) * ch + k))); \

    a_ptr += a->step; \

    b_ptr += db->step; \

    for (i = 1; i < a->rows - 1; 
i++498
) \

    { \

      for (k = 0; k < ch; 
k++498
) \

        _for_set(b_ptr, k, 2 * (_for_get(a_ptr + a->step, ch + k) - _for_get(a_ptr, k))); \

      for (j = 1; j < a->cols - 1; 
j++248k
) \

        
for (k = 0; 248k
k < ch; 
k++248k
) \

          
_for_set498
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, (j + 1) * ch + k) - _for_get(a_ptr - a->step, (j - 1) * ch + k)); \

      for (k = 0; k < ch; 
k++498
) \

        _for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k) - _for_get(a_ptr - a->step, (a->cols - 2) * ch + k))); \

      a_ptr += a->step; \

      b_ptr += db->step; \

    } \

    for (k = 0; k < ch; 
k++1
) \

      _for_set(b_ptr, k, 2 * (_for_get(a_ptr, k) - _for_get(a_ptr - a->step, k))); \

    for (j = 1; j < a->cols; 
j++499
) \

      
for (k = 0; 499
k < ch; 
k++499
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k) - _for_get(a_ptr - a->step, (j - 1) * ch + k)));

    ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);

#undef for_block

  } else if ((dx == 1 && 
dy == -10
) || (dx == -1 && 
dy == 11
)) {

    /* 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; 1
k < ch; 
k++1
) \

      _for_set(b_ptr, k, 2 * (_for_get(a_ptr + a->step, k) - _for_get(a_ptr, k))); \

    for (j = 1; j < a->cols; 
j++499
) \

      
for (k = 0; 499
k < ch; 
k++499
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr + a->step, (j - 1) * ch + k) - _for_get(a_ptr, j * ch + k))); \

    a_ptr += a->step; \

    b_ptr += db->step; \

    for (i = 1; i < a->rows - 1; 
i++498
) \

    { \

      for (k = 0; k < ch; 
k++498
) \

        _for_set(b_ptr, k, 2 * (_for_get(a_ptr, k) - _for_get(a_ptr - a->step, ch + k))); \

      for (j = 1; j < a->cols - 1; 
j++248k
) \

        
for (k = 0; 248k
k < ch; 
k++248k
) \

          
_for_set498
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, (j - 1) * ch + k) - _for_get(a_ptr - a->step, (j + 1) * ch + k)); \

      for (k = 0; k < ch; 
k++498
) \

        _for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr + a->step, (a->cols - 2) * ch + k) - _for_get(a_ptr, (a->cols - 1) * ch + k))); \

      a_ptr += a->step; \

      b_ptr += db->step; \

    } \

    for (j = 0; j < a->cols - 1; 
j++499
) \

      
for (k = 0; 499
k < ch; 
k++499
) \

        
_for_set1
(b_ptr, j * ch + k, 2 * (_for_get(a_ptr, j * ch + k) - _for_get(a_ptr - a->step, (j + 1) * ch + k))); \

    for (k = 0; k < ch; 
k++1
) \

      _for_set(b_ptr, (a->cols - 1) * ch + k, 2 * (_for_get(a_ptr, (a->cols - 1) * ch + k) - _for_get(a_ptr - a->step, (a->cols - 1) * ch + k)));

    ccv_matrix_getter(a->type, ccv_matrix_setter, db->type, for_block);

#undef for_block

  } else if (dx == 3 && 
dy == 02
) {

    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; 2
j < a->cols; 
j++1.00k
) \

      
for (k = 0; 1.00k
k < ch; 
k++1.00k
) \

        
_for_set_b2
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k) + 3 * _for_get(a_ptr, j * ch + k)); \

    a_ptr += a->step; \

    b_ptr += db->step; \

    for (i = 1; i < a->rows - 1; 
i++996
) \

    { \

      for (j = 0; j < a->cols; 
j++498k
) \

        
for (k = 0; 498k
k < ch; 
k++498k
) \

          
_for_set_b996
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k) + 2 * _for_get(a_ptr, j * ch + k) + _for_get(a_ptr - a->step, j * ch + k)); \

      a_ptr += a->step; \

      b_ptr += db->step; \

    } \

    for (j = 0; j < a->cols; 
j++1.00k
) \

      
for (k = 0; 1.00k
k < ch; 
k++1.00k
) \

        
_for_set_b2
(b_ptr, j * ch + k, 3 * _for_get(a_ptr, j * ch + k) + _for_get(a_ptr - a->step, j * ch + k)); \

    b_ptr = db->data.u8; \

    for (i = 0; i < a->rows; 
i++1.00k
) \

    { \

      for (k = 0; k < ch; 
k++1.00k
) \

        _for_set_b(buf, k, _for_get_b(b_ptr, ch + k) - _for_get_b(b_ptr, k)); \

      for (j = 1; j < a->cols - 1; 
j++498k
) \

        
for (k = 0; 498k
k < ch; 
k++498k
) \

          
_for_set_b1.00k
(buf, j * ch + k, _for_get_b(b_ptr, (j + 1) * ch + k) - _for_get_b(b_ptr, (j - 1) * ch + k)); \

      for (k = 0; k < ch; 
k++1.00k
) \

        _for_set_b(buf, (a->cols - 1) * ch + k, _for_get_b(b_ptr, (a->cols - 1) * ch + k) - _for_get_b(b_ptr, (a->cols - 2) * ch + k)); \

      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 == 33
) {

    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; 2
j < a->cols; 
j++1.00k
) \

      
for (k = 0; 1.00k
k < ch; 
k++1.00k
) \

        
_for_set_b2
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k) - _for_get(a_ptr, j * ch + k)); \

    a_ptr += a->step; \

    b_ptr += db->step; \

    for (i = 1; i < a->rows - 1; 
i++996
) \

    { \

      for (j = 0; j < a->cols; 
j++498k
) \

        
for (k = 0; 498k
k < ch; 
k++498k
) \

          
_for_set_b996
(b_ptr, j * ch + k, _for_get(a_ptr + a->step, j * ch + k) - _for_get(a_ptr - a->step, j * ch + k)); \

      a_ptr += a->step; \

      b_ptr += db->step; \

    } \

    for (j = 0; j < a->cols; 
j++1.00k
) \

      
for (k = 0; 1.00k
k < ch; 
k++1.00k
) \

        
_for_set_b2
(b_ptr, j * ch + k, _for_get(a_ptr, j * ch + k) - _for_get(a_ptr - a->step, j * ch + k)); \

    b_ptr = db->data.u8; \

    for (i = 0; i < a->rows; 
i++1.00k
) \

    { \

      for (k = 0; k < ch; 
k++1.00k
) \

        _for_set_b(buf, k, _for_get_b(b_ptr, ch + k) + 3 * _for_get_b(b_ptr, k)); \

      for (j = 1; j < a->cols - 1; 
j++498k
) \

        
for (k = 0; 498k
k < ch; 
k++498k
) \

          
_for_set_b1.00k
(buf, j * ch + k, _for_get_b(b_ptr, (j + 1) * ch + k) + 2 * _for_get_b(b_ptr, j * ch + k) + _for_get_b(b_ptr, (j - 1) * ch + k)); \

      for (k = 0; k < ch; 
k++1.00k
) \

        _for_set_b(buf, (a->cols - 1) * ch + k, _for_get_b(b_ptr, (a->cols - 2) * ch + k) + 3 * _for_get_b(b_ptr, (a->cols - 1) * ch + k)); \

      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++10
)

    {

      ((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++10
)

      {

        // 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; 2
i < a->rows; 
i++1.00k
) \

    { \

      for (j = 0; j < hfz; 
j++2.00k
) \

        
for (k = 0; 2.00k
k < ch; 
k++2.00k
) \

          
_for_set_b1.00k
(buf, j * ch + k, _for_get(a_ptr, k)); \

      for (j = 0; j < a->cols; 
j++500k
) \

        
for (k = 0; 500k
k < ch; 
k++500k
) \

          
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, j * ch + k)); \

      for (j = a->cols; j < a->cols + hfz; 
j++2.00k
) \

        
for (k = 0; 2.00k
k < ch; 
k++2.00k
) \

          
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, (a->cols - 1) * ch + k)); \

      for (j = 0; j < a->cols; 
j++500k
) \

      { \

        for (c = 0; c < ch; 
c++500k
) \

        { \

          _for_type_b sum = 0; \

          for (k = 0; k < fsz; 
k++2.50M
) \

            sum += _for_get_b(buf, (j + k) * ch + c) * 
_for_get_b500k
(df, k); \

          _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++1.00k
) \

    { \

      for (j = 0; j < hfz; 
j++2.00k
) \

        
for (k = 0; 2.00k
k < ch; 
k++2.00k
) \

          
_for_set_b1.00k
(buf, j * ch + k, _for_get_b(b_ptr, i * ch + k)); \

      for (j = 0; j < a->rows; 
j++500k
) \

        
for (k = 0; 500k
k < ch; 
k++500k
) \

          
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + j * db->step, i * ch + k)); \

      for (j = a->rows; j < a->rows + hfz; 
j++2.00k
) \

        
for (k = 0; 2.00k
k < ch; 
k++2.00k
) \

          
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + (a->rows - 1) * db->step, i * ch + k)); \

      for (j = 0; j < a->rows; 
j++500k
) \

      { \

        for (c = 0; c < ch; 
c++500k
) \

        { \

          _for_type_b sum = 0; \

          for (k = 0; k < fsz; 
k++2.50M
) \

            sum += _for_get_b(buf, (j + k) * ch + c) * 
_for_get_b500k
(gf, k); \

          _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++500
)

  {

    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++2.40M
)

  {

    for (j = 0; j < a->cols / 2; 
j++38.7M
)

    {

      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) ? 
00
 : ccv_cache_generate_signature("ccv_flip_y", 10, sig, CCV_EOF_SIGN);

  if (type & CCV_FLIP_X)

    sig = (a->sig == 0) ? 
075.1k
 : 
ccv_cache_generate_signature("ccv_flip_x", 10, sig, 2
CCV_EOF_SIGN2
);

  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_TYPE0
(type) | 0
CCV_GET_CHANNEL0
(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; 1
i < fsz; 
i++27
)

    tw += ((double*)filter)[i] = exp(-((i - hfz) * (i - hfz)) / (2.0 * sigma * sigma));

  int no_8u_type = (db->type & CCV_8U) ? CCV_32S : 
db->type0
;

  if (no_8u_type & CCV_32S)

  {

    tw = 256.0 / tw;

    for (i = 0; i < fsz; 
i++27
)

      ((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; 1
i < a->rows; 
i++600
) \

  { \

    for (j = 0; j < hfz; 
j++7.80k
) \

      
for (k = 0; 7.80k
k < ch; 
k++23.4k
) \

        
_for_set_b600
(buf, j * ch + k, _for_get_a(a_ptr, k)); \

    for (j = 0; j < a->cols * ch; 
j++1.44M
) \

      _for_set_b(buf, j + hfz * ch, _for_get_a(a_ptr, j)); \

    for (j = a->cols; j < hfz + a->cols; 
j++7.80k
) \

      
for (k = 0; 7.80k
k < ch; 
k++23.4k
) \

        
_for_set_b600
(buf, j * ch + hfz * ch + k, _for_get_a(a_ptr, (a->cols - 1) * ch + k)); \

    for (j = 0; j < a->cols * ch; 
j++1.44M
) \

    { \

      _for_type sum = 0; \

      for (k = 0; k < fsz; 
k++38.8M
) \

        sum += _for_get_b(buf, k * ch + j) * 
_for_get_b1.44M
(filter, k); \

      _for_set_b(buf, j, sum, 8); \

    } \

    for (j = 0; j < a->cols * ch; 
j++1.44M
) \

      _for_set_a(b_ptr, j, _for_get_b(buf, j)); \

    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; 1
i < a->cols * ch; 
i++2.40k
) \

  { \

    for (j = 0; j < hfz; 
j++31.2k
) \

      _for_set_b(buf, j, _for_get_a(b_ptr, i)); \

    for (j = 0; j < a->rows; 
j++1.44M
) \

      _for_set_b(buf, j + hfz, _for_get_a(b_ptr + j * db->step, i)); \

    for (j = a->rows; j < hfz + a->rows; 
j++31.2k
) \

      _for_set_b(buf, j + hfz, _for_get_a(b_ptr + (a->rows - 1) * db->step, i)); \

    for (j = 0; j < a->rows; 
j++1.44M
) \

    { \

      _for_type sum = 0; \

      for (k = 0; k < fsz; 
k++38.8M
) \

        sum += _for_get_b(buf, k + j) * 
_for_get_b1.44M
(filter, k); \

      _for_set_b(buf, j, sum, 8); \

    } \

    for (j = 0; j < a->rows; 
j++1.44M
) \

      _for_set_a(b_ptr + j * db->step, i, _for_get_b(buf, j)); \

  }

  ccv_matrix_typeof_setter_getter(no_8u_type, ccv_matrix_setter_getter, db->type, for_block);

#undef for_block

}

void ccv_erode(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int fsz)

{

  ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_erode(%d)", fsz), a->sig, CCV_EOF_SIGN);

  type = (type == 0) ? CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type) : 
CCV_GET_DATA_TYPE0
(type) | 0
CCV_GET_CHANNEL0
(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* buf = (unsigned char*)alloca(sizeof(double) * ch * (fsz + ccv_max(a->rows, a->cols)));

  const int hfz = fsz / 2;

  unsigned char* a_ptr = a->data.u8;

  unsigned char* b_ptr = db->data.u8;

#define for_block(_for_get, _for_type_b, _for_set_b, _for_get_b) \

  
for (i = 0; 2
i < a->rows; 
i++1.00k
) \

  { \

    for (j = 0; j < hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, j * ch + k, _for_get(a_ptr, k)); \

    for (j = 0; j < a->cols; 
j++500k
) \

      
for (k = 0; 500k
k < ch; 
k++1.50M
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, j * ch + k)); \

    for (j = a->cols; j < a->cols + hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, (a->cols - 1) * ch + k)); \

    for (j = 0; j < a->cols; 
j++500k
) \

    { \

      for (c = 0; c < ch; 
c++1.50M
) \

      { \

        _for_type_b minimum = _for_get_b(buf, j * ch + c); \

        for (k = 1; k < fsz; 
k++4.50M
) \

          minimum = ccv_min(minimum, _for_get_b(buf, (j + k) * ch + c)); \

        _for_set_b(b_ptr, j * ch + c, minimum); \

      } \

    } \

    a_ptr += a->step; \

    b_ptr += db->step; \

  } \

  b_ptr = db->data.u8; \

  for (i = 0; i < a->cols; 
i++1.00k
) \

  { \

    for (j = 0; j < hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, j * ch + k, _for_get_b(b_ptr, i * ch + k)); \

    for (j = 0; j < a->rows; 
j++500k
) \

      
for (k = 0; 500k
k < ch; 
k++1.50M
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + j * db->step, i * ch + k)); \

    for (j = a->rows; j < a->rows + hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + (a->rows - 1) * db->step, i * ch + k)); \

    for (j = 0; j < a->rows; 
j++500k
) \

    { \

      for (c = 0; c < ch; 
c++1.50M
) \

      { \

        _for_type_b minimum = _for_get_b(buf, j * ch + c); \

        for (k = 1; k < fsz; 
k++4.50M
) \

          minimum = ccv_min(minimum, _for_get_b(buf, (j + k) * ch + c)); \

        _for_set_b(b_ptr + j * db->step, i * ch + c, minimum); \

      } \

    } \

  }

  ccv_matrix_getter(a->type, ccv_matrix_typeof_setter_getter, db->type, for_block);

#undef for_block

}

void ccv_dilate(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int fsz)

{

  ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_dilate(%d)", fsz), a->sig, CCV_EOF_SIGN);

  type = (type == 0) ? CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type) : 
CCV_GET_DATA_TYPE0
(type) | 0
CCV_GET_CHANNEL0
(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* buf = (unsigned char*)alloca(sizeof(double) * ch * (fsz + ccv_max(a->rows, a->cols)));

  const int hfz = fsz / 2;

  unsigned char* a_ptr = a->data.u8;

  unsigned char* b_ptr = db->data.u8;

#define for_block(_for_get, _for_type_b, _for_set_b, _for_get_b) \

  
for (i = 0; 2
i < a->rows; 
i++1.00k
) \

  { \

    for (j = 0; j < hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, j * ch + k, _for_get(a_ptr, k)); \

    for (j = 0; j < a->cols; 
j++500k
) \

      
for (k = 0; 500k
k < ch; 
k++1.50M
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, j * ch + k)); \

    for (j = a->cols; j < a->cols + hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get(a_ptr, (a->cols - 1) * ch + k)); \

    for (j = 0; j < a->cols; 
j++500k
) \

    { \

      for (c = 0; c < ch; 
c++1.50M
) \

      { \

        _for_type_b maximum = _for_get_b(buf, j * ch + c); \

        for (k = 1; k < fsz; 
k++4.50M
) \

          maximum = ccv_max(maximum, _for_get_b(buf, (j + k) * ch + c)); \

        _for_set_b(b_ptr, j * ch + c, maximum); \

      } \

    } \

    a_ptr += a->step; \

    b_ptr += db->step; \

  } \

  b_ptr = db->data.u8; \

  for (i = 0; i < a->cols; 
i++1.00k
) \

  { \

    for (j = 0; j < hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, j * ch + k, _for_get_b(b_ptr, i * ch + k)); \

    for (j = 0; j < a->rows; 
j++500k
) \

      
for (k = 0; 500k
k < ch; 
k++1.50M
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + j * db->step, i * ch + k)); \

    for (j = a->rows; j < a->rows + hfz; 
j++1.50k
) \

      
for (k = 0; 1.50k
k < ch; 
k++4.50k
) \

        
_for_set_b1.00k
(buf, (j + hfz) * ch + k, _for_get_b(b_ptr + (a->rows - 1) * db->step, i * ch + k)); \

    for (j = 0; j < a->rows; 
j++500k
) \

    { \

      for (c = 0; c < ch; 
c++1.50M
) \

      { \

        _for_type_b maximum = _for_get_b(buf, j * ch + c); \

        for (k = 1; k < fsz; 
k++4.50M
) \

          maximum = ccv_max(maximum, _for_get_b(buf, (j + k) * ch + c)); \

        _for_set_b(b_ptr + j * db->step, i * ch + c, maximum); \

      } \

    } \

  }

  ccv_matrix_getter(a->type, ccv_matrix_typeof_setter_getter, db->type, for_block);

#undef for_block

}