Skip to content

Commit

Permalink
added ccv_distance_transform which implements generalized squared
Browse files Browse the repository at this point in the history 
euclidean distance transform in linear time.

added test cases for ccv_convolve and ccv_distane_transform
  • Loading branch information
@liuliu
liuliu committed Mar 29, 2012
1 parent eb1483b commit 5fe7f83
Show file tree
Hide file tree
Showing 15 changed files with 450 additions and 266 deletions.
2 changes: 0 additions & 2 deletions bin/dpmdetect.c
View file
Expand Up @@ -22,7 +22,6 @@ int main(int argc, char** argv)
unsigned int elapsed_time = get_current_time();
ccv_dpm_param_t params = { .interval = 5, .min_neighbors = 2, .flags = 0, .size = ccv_size(root_classifier->root.size.width * 8, root_classifier->root.size.height * 8) };
ccv_array_t* seq = ccv_dpm_detect_objects(image, &root_classifier, 1, params);
/*
elapsed_time = get_current_time() - elapsed_time;
for (i = 0; i < seq->rnum; i++)
{
Expand All @@ -31,7 +30,6 @@ int main(int argc, char** argv)
}
printf("total : %d in time %dms\n", seq->rnum, elapsed_time);
ccv_array_free(seq);
*/
ccv_matrix_free(image);
}
ccv_drain_cache();
Expand Down
2 changes: 1 addition & 1 deletion bin/makefile
View file
Expand Up @@ -2,7 +2,7 @@ CC = clang
LDFLAGS = -L"../lib" -lccv -pthread -ljpeg -lpng -lfftw3 -lz -lgsl -lblas -lm # -lgomp # -lOpenCL -L"/opt/ati-stream-sdk/lib/x86_64"
CXXFLAGS = -O3 -msse2 -Wall -I"../lib" # -fopenmp -D USE_OPENMP # -D USE_OPENCL

TARGETS = bbffmt siftmatch bbfcreate bbfdetect swtdetect swtcreate convert distance blur
TARGETS = bbffmt siftmatch bbfcreate bbfdetect swtdetect swtcreate convert

all: libccv.a $(TARGETS)

Expand Down
30 changes: 14 additions & 16 deletions lib/ccv.h
View file
Expand Up @@ -14,6 +14,7 @@
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <xmmintrin.h>
#include <assert.h>
Expand Down Expand Up @@ -81,6 +82,7 @@ typedef struct {
float fl;
int64_t l;
double db;
void* ptr;
} tag;
ccv_matrix_cell_t data;
} ccv_dense_matrix_t;
Expand Down Expand Up @@ -337,15 +339,7 @@ void ccv_enable_cache(size_t size);
case CCV_64F: ((double*)(ptr))[(i)] = (double)value; break; \
default: ((unsigned char*)(ptr))[(i)] = ccv_clamp((int)(value) >> factor, 0, 255); }


/* unswitch for loop macros */
/* the new added macro in order to do for loop expansion in a way that, you can
* expand a for loop by inserting different code snippet */
#define ccv_unswitch_block(param, block, ...) { block(__VA_ARGS__, param); }
#define ccv_unswitch_block_a(param, block, ...) { block(__VA_ARGS__, param); }
#define ccv_unswitch_block_b(param, block, ...) { block(__VA_ARGS__, param); }
/* the factor used to provide higher accuracy in integer type (all integer
* computation in some cases) */
/* the factor used to provide higher accuracy in integer type (all integer computation in some cases) */
#define _ccv_get_32s_value(ptr, i, factor) (((int*)(ptr))[(i)] << factor)
#define _ccv_get_32f_value(ptr, i, factor) ((float*)(ptr))[(i)]
#define _ccv_get_64s_value(ptr, i, factor) (((int64_t*)(ptr))[(i)] << factor)
Expand Down Expand Up @@ -533,8 +527,11 @@ int ccv_write(ccv_dense_matrix_t* mat, char* out, int* len, int type, void* conf
double ccv_trace(ccv_matrix_t* mat);

enum {
CCV_L2_NORM = 0x01,
CCV_L1_NORM = 0x02,
CCV_L2_NORM = 0x01, // |dx| + |dy|
CCV_L1_NORM = 0x02, // sqrt(dx^2 + dy^2)
CCV_GSEDT = 0x04, // Generalized Squared Euclidean Distance Transform:
// a * dx + b * dy + c * dx^2 + d * dy^2, when combined with CCV_L1_NORM:
// a * |dx| + b * |dy| + c * dx^2 + d * dy^2
};

enum {
Expand Down Expand Up @@ -676,20 +673,21 @@ typedef double(*ccv_convolve_kernel_f)(double x, double y, void*);
void ccv_convolve_kernel(ccv_dense_matrix_t* x, ccv_convolve_kernel_f func, void* data);

/* modern numerical algorithms */
void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, double dx, double dy, double dxx, double dyy, int flag);
void ccv_sparse_coding(ccv_matrix_t* x, int k, ccv_matrix_t** A, int typeA, ccv_matrix_t** y, int typey);
void ccv_compressive_sensing_reconstruct(ccv_matrix_t* a, ccv_matrix_t* x, ccv_matrix_t** y, int type);

/* basic computer vision algorithms / or build blocks */
void ccv_sobel(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int dx, int dy);
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);
void ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int size);
void ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int b_type, int sbin, int size);
void ccv_canny(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int size, double low_thresh, double high_thresh);

enum {
CCV_INTER_AREA = 0x01,
CCV_INTER_LINEAR = 0X02,
CCV_INTER_CUBIC = 0X03,
CCV_INTER_LACZOS = 0X04,
CCV_INTER_AREA = 0x01,
CCV_INTER_LINEAR = 0X02,
CCV_INTER_CUBIC = 0X03,
CCV_INTER_LANCZOS = 0X04,
};

void ccv_resample(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int btype, int rows, int cols, int type);
Expand Down
198 changes: 154 additions & 44 deletions lib/ccv_basic.c
View file
Expand Up @@ -244,7 +244,7 @@ static void _ccv_atan2(float* x, float* y, float* angle, float* mag, int len)

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)
{
assert(a->type & CCV_C1);
assert(CCV_GET_CHANNEL(a->type) == CCV_C1);
ccv_declare_matrix_signature(tsig, a->sig != 0, ccv_sign_with_literal("ccv_gradient_theta"), a->sig, 0);
ccv_declare_matrix_signature(msig, a->sig != 0, ccv_sign_with_literal("ccv_gradient_m"), a->sig, 0);
ccv_dense_matrix_t* dtheta = *theta = ccv_dense_matrix_renew(*theta, a->rows, a->cols, CCV_32F | CCV_C1, CCV_32F | CCV_C1, tsig);
Expand All @@ -260,53 +260,164 @@ void ccv_gradient(ccv_dense_matrix_t* a, ccv_dense_matrix_t** theta, int ttype,
ccv_matrix_free(ty);
}

void ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int size)
void ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int b_type, int sbin, int size)
{
assert(a->type & CCV_C1);
int border_size = size / 2;
ccv_declare_matrix_signature(sig, a->sig != 0, ccv_sign_with_literal("ccv_hog"), a->sig, 0);
type = (type == 0) ? CCV_32S | CCV_C1 : CCV_GET_DATA_TYPE(type) | CCV_C1;
ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, a->rows - border_size * 2, (a->cols - border_size * 2) * 8, CCV_C1 | CCV_ALL_DATA_TYPE, type, sig);
assert(CCV_GET_CHANNEL(a->type) == CCV_C1);
int rows = a->rows / size;
int cols = a->cols / size;
ccv_declare_matrix_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_hog(%d,%d)", sbin, size), a->sig, 0);
ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, rows, cols, CCV_64F | (4 + sbin * 3), CCV_64F | (4 + sbin * 3), sig);
ccv_matrix_return_if_cached(, db);
ccv_dense_matrix_t* ag = 0;
ccv_dense_matrix_t* mg = 0;
ccv_gradient(a, &ag, 0, &mg, 0, 3, 3);
int i, j, x, y;
ag->type = (ag->type & ~CCV_32F) | CCV_32S;
for (i = 0; i < a->rows * a->cols; i++)
ag->data.i[i] = ((int)(ag->data.fl[i] / 45 + 0.5)) & 0x7;
int* agi = ag->data.i;
float* mgfl = mg->data.fl;
int* bi = db->data.i;
for (y = 0; y <= a->rows - size; y++)
ccv_gradient(a, &ag, 0, &mg, 0, 1, 1);
float* agp = ag->data.fl;
float* mgp = mg->data.fl;
int i, j, k;
ccv_dense_matrix_t* cn = ccv_dense_matrix_new(rows, cols, CCV_64F | (sbin * 2), 0, 0);
ccv_dense_matrix_t* ca = ccv_dense_matrix_new(rows, cols, CCV_64F | CCV_C1, 0, 0);
ccv_zero(cn);
double* cnp = cn->data.db;
int sizec = 0;
for (i = 0; i < rows * size; i++)
{
float hog[] = { 0, 0, 0, 0, 0, 0, 0, 0 };
for (i = 0; i < size; i++)
for (j = 0; j < size; j++)
hog[agi[i * a->cols + j]] += mgfl[i * a->cols + j];
for (i = 0; i < 8; i++)
bi[i] = (int)hog[i];
bi += 8;
mgfl++;
agi++;
for (x = 1; x <= a->cols - size; x++)
for (j = 0; j < cols; j++)
{
for (i = 0; i < size; i++)
for (k = j * size; k < j * size + size; k++)
{
hog[agi[i * a->cols - 1]] -= mgfl[i * a->cols - 1];
hog[agi[i * a->cols - 1 + size]] += mgfl[i * a->cols - 1 + size];
int ag0, ag1;
double agr;
agr = (agp[k] / 360.0) * (sbin * 2);
ag0 = (int)agr;
ag1 = (ag0 + 1 < sbin * 2) ? ag0 + 1 : 0;
agr = agr - ag0;
cnp[ag0] += (1.0 - agr) * mgp[k] / 255.0;
cnp[ag1] += agr * mgp[k] / 255.0;
}
for (i = 0; i < 8; i++)
bi[i] = (int)hog[i];
bi += 8;
mgfl++;
agi++;
cnp += 2 * sbin;
}
agi += border_size * 2;
mgfl += border_size * 2;
agp += a->cols;
mgp += a->cols;
if (++sizec < size)
cnp -= cn->cols;
else
sizec = 0;
}
ccv_matrix_free(ag);
ccv_matrix_free(mg);
cnp = cn->data.db;
double* cap = ca->data.db;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
*cap = 0;
for (k = 0; k < sbin * 2; k++)
{
*cap += (*cnp) * (*cnp);
cnp++;
}
cap++;
}
}
cnp = cn->data.db;
cap = ca->data.db;
ccv_zero(db);
double* dbp = db->data.db;
// normalize sbin direction-sensitive and sbin * 2 insensitive over 4 normalization factor
// accumulating them over sbin * 2 + sbin + 4 channels
// TNA - truncation - normalization - accumulation
#define TNA(idx, a, b, c, d) \
{ \
double norm = 1.0 / sqrt(cap[a] + cap[b] + cap[c] + cap[d] + 1e-4); \
for (k = 0; k < sbin * 2; k++) \
{ \
double v = 0.5 * ccv_min(cnp[k] * norm, 0.2); \
dbp[5 + sbin + k] += v; \
dbp[1 + idx] += v; \
} \
dbp[sbin * 3 + idx] *= 0.2357; \
for (k = 0; k < sbin; k++) \
{ \
double v = 0.5 * ccv_min((cnp[k] + cnp[k + sbin]) * norm, 0.2); \
dbp[5 + k] += v; \
} \
}
TNA(0, 0, 0, 0, 0);
TNA(1, 1, 1, 0, 0);
TNA(2, 0, ca->cols, ca->cols, 0);
TNA(3, 1, ca->cols + 1, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
for (j = 1; j < cols - 1; j++)
{
TNA(0, -1, -1, 0, 0);
TNA(1, 1, 1, 0, 0);
TNA(2, -1, ca->cols - 1, ca->cols, 0);
TNA(3, 1, ca->cols + 1, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
}
TNA(0, -1, -1, 0, 0);
TNA(1, 0, 0, 0, 0);
TNA(2, -1, ca->cols - 1, ca->cols, 0);
TNA(3, 0, ca->cols, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
for (i = 1; i < rows - 1; i++)
{
TNA(0, 0, -ca->cols, -ca->cols, 0);
TNA(1, 1, -ca->cols + 1, -ca->cols, 0);
TNA(2, 0, ca->cols, ca->cols, 0);
TNA(3, 1, ca->cols + 1, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
for (j = 1; j < cols - 1; j++)
{
TNA(0, -1, -ca->cols - 1, -ca->cols, 0);
TNA(1, 1, -ca->cols + 1, -ca->cols, 0);
TNA(2, -1, ca->cols - 1, ca->cols, 0);
TNA(3, 1, ca->cols + 1, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
}
TNA(0, -1, -ca->cols - 1, -ca->cols, 0);
TNA(1, 0, -ca->cols, -ca->cols, 0);
TNA(2, -1, ca->cols - 1, ca->cols, 0);
TNA(3, 0, ca->cols, ca->cols, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
}
TNA(0, 0, -ca->cols, -ca->cols, 0);
TNA(1, 1, -ca->cols + 1, -ca->cols, 0);
TNA(2, 0, 0, 0, 0);
TNA(3, 1, 1, 0, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
for (j = 1; j < cols - 1; j++)
{
TNA(0, -1, -ca->cols - 1, -ca->cols, 0);
TNA(1, 1, -ca->cols + 1, -ca->cols, 0);
TNA(2, -1, -1, 0, 0);
TNA(3, 1, 1, 0, 0);
cnp += 2 * sbin;
dbp += 3 * sbin + 4;
cap++;
}
TNA(0, -1, -ca->cols - 1, -ca->cols, 0);
TNA(1, 0, -ca->cols, -ca->cols, 0);
TNA(2, -1, -1, 0, 0);
TNA(3, 0, 0, 0, 0);
#undef TNA
ccv_matrix_free(cn);
ccv_matrix_free(ca);
}

/* it is a supposely cleaner and faster implementation than original OpenCV (ccv_canny_deprecated,
Expand Down Expand Up @@ -678,7 +789,7 @@ void ccv_resample(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int btype, int
break;
case CCV_INTER_CUBIC:
break;
case CCV_INTER_LACZOS:
case CCV_INTER_LANCZOS:
break;
}
}
Expand Down Expand Up @@ -707,7 +818,7 @@ void ccv_sample_down(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, in
* it is not desirable to have that offset when we try to wrap it into our 5-row buffer (
* because in later rearrangement, we have no src_y to backup the arrangement). In
* such micro scope, we managed to stripe 5 addition into one shift and addition. */
#define for_block(x_block, _for_get_a, _for_get, _for_set, _for_set_b) \
#define for_block(_for_get_a, _for_get, _for_set, _for_set_b) \
for (dy = 0; dy < db->rows; dy++) \
{ \
for(; sy <= dy * 2 + 2 + src_y; sy++) \
Expand All @@ -730,20 +841,19 @@ void ccv_sample_down(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, in
b_ptr += db->step; \
}
int no_8u_type = (a->type & CCV_8U) ? CCV_32S : a->type;
/* here is the new technique to expand for loop with condition in manual way */
if (src_x > 0)
{
#define x_block(_for_get_a, _for_get, _for_set, _for_set_b) \
for (dx = cols0 * ch; dx < db->cols * ch; dx += ch) \
for (k = 0; k < ch; k++) \
_for_set(row, dx + k, _for_get_a(a_ptr, tab[dx * 2 + sx + k], 0) * 6 + (_for_get_a(a_ptr, tab[dx * 2 + sx + k - ch], 0) + _for_get_a(a_ptr, tab[dx * 2 + sx + k + ch], 0)) * 4 + _for_get_a(a_ptr, tab[dx * 2 + sx + k - ch * 2], 0) + _for_get_a(a_ptr, tab[dx * 2 + sx + k + ch * 2], 0), 0);
ccv_unswitch_block(x_block, ccv_matrix_getter_a, a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
ccv_matrix_getter_a(a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
#undef x_block
} else {
#define x_block(_for_get_a, _for_get, _for_set, _for_set_b) \
for (k = 0; k < ch; k++) \
_for_set(row, (db->cols - 1) * ch + k, _for_get_a(a_ptr, a->cols * ch + sx - ch + k, 0) * 10 + _for_get_a(a_ptr, (a->cols - 2) * ch + sx + k, 0) * 5 + _for_get_a(a_ptr, (a->cols - 3) * ch + sx + k, 0), 0);
ccv_unswitch_block(x_block, ccv_matrix_getter_a, a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
ccv_matrix_getter_a(a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
#undef x_block
}
#undef for_block
Expand All @@ -767,7 +877,7 @@ void ccv_sample_up(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int
int bufstep = db->cols * ch * ccv_max(CCV_GET_DATA_TYPE_SIZE(db->type), sizeof(int));
unsigned char* b_ptr = db->data.ptr;
/* why src_y * 2: the same argument as in ccv_sample_down */
#define for_block(x_block, _for_get_a, _for_get, _for_set, _for_set_b) \
#define for_block(_for_get_a, _for_get, _for_set, _for_set_b) \
for (y = 0; y < a->rows; y++) \
{ \
for (; sy <= y + 1 + src_y; sy++) \
Expand Down Expand Up @@ -820,7 +930,7 @@ void ccv_sample_up(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int
_for_set(row, x * 2 + k, _for_get_a(a_ptr, tab[x + sx - ch + k], 0) + _for_get_a(a_ptr, tab[x + sx + k], 0) * 6 + _for_get_a(a_ptr, tab[x + sx + ch + k], 0), 0); \
_for_set(row, x * 2 + ch + k, (_for_get_a(a_ptr, tab[x + sx + k], 0) + _for_get_a(a_ptr, tab[x + sx + ch + k], 0)) * 4, 0); \
}
ccv_unswitch_block(x_block, ccv_matrix_getter_a, a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
ccv_matrix_getter_a(a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
#undef x_block
} else {
#define x_block(_for_get_a, _for_get, _for_set, _for_set_b) \
Expand All @@ -829,7 +939,7 @@ void ccv_sample_up(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int
_for_set(row, (a->cols - 1) * 2 * ch + k, _for_get_a(a_ptr, (a->cols - 2) * ch + k, 0) + _for_get_a(a_ptr, (a->cols - 1) * ch + k, 0) * 7, 0); \
_for_set(row, (a->cols - 1) * 2 * ch + ch + k, _for_get_a(a_ptr, (a->cols - 1) * ch + k, 0) * 4, 0); \
}
ccv_unswitch_block(x_block, ccv_matrix_getter_a, a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
ccv_matrix_getter_a(a->type, ccv_matrix_getter, no_8u_type, ccv_matrix_setter, no_8u_type, ccv_matrix_setter_b, db->type, for_block);
#undef x_block
}
#undef for_block
Expand Down

0 comments on commit 5fe7f83

Please sign in to comment.