diff --git a/bin/dpmdetect.c b/bin/dpmdetect.c index f02c15441..d3b88ca2e 100644 --- a/bin/dpmdetect.c +++ b/bin/dpmdetect.c @@ -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++) { @@ -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(); diff --git a/bin/makefile b/bin/makefile index d4733a2c1..c42ea112a 100644 --- a/bin/makefile +++ b/bin/makefile @@ -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) diff --git a/lib/ccv.h b/lib/ccv.h index 4b7ca1aa6..7c76aa654 100644 --- a/lib/ccv.h +++ b/lib/ccv.h @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -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; @@ -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) @@ -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 { @@ -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); diff --git a/lib/ccv_basic.c b/lib/ccv_basic.c index 12dc9dac6..c9dbe1329 100644 --- a/lib/ccv_basic.c +++ b/lib/ccv_basic.c @@ -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); @@ -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, @@ -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; } } @@ -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++) \ @@ -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 @@ -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++) \ @@ -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) \ @@ -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 diff --git a/lib/ccv_dpm.c b/lib/ccv_dpm.c index dd9d68299..50bdee7ac 100644 --- a/lib/ccv_dpm.c +++ b/lib/ccv_dpm.c @@ -4,167 +4,6 @@ void ccv_dpm_classifier_lsvm_new(ccv_dense_matrix_t** posimgs, int posnum, char* { } -static inline float _ccv_lsvm_match(ccv_dense_matrix_t* a, ccv_dense_matrix_t* tpl, ccv_dense_matrix_t** b) -{ -} - -// this is specific HOG computation for dpm, I may later change it to ccv_hog -static void _ccv_hog(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int b_type, int sbin, int size) -{ - ccv_dense_matrix_t* ag = 0; - ccv_dense_matrix_t* mg = 0; - ccv_gradient(a, &ag, 0, &mg, 0, 1, 1); - float* agp = ag->data.fl; - float* mgp = mg->data.fl; - int i, j, k; - int rows = a->rows / size; - int cols = a->cols / size; - ccv_dense_matrix_t* cn = ccv_dense_matrix_new(rows, cols, CCV_32F | (sbin * 2), 0, 0); - ccv_dense_matrix_t* ca = ccv_dense_matrix_new(rows, cols, CCV_64F | CCV_C1, 0, 0); - ccv_dense_matrix_t* db = *b = ccv_dense_matrix_new(rows, cols, CCV_32F | (4 + sbin * 3), 0, 0); - ccv_zero(cn); - float* cnp = cn->data.fl; - int sizec = 0; - for (i = 0; i < rows * size; i++) - { - for (j = 0; j < cols; j++) - { - for (k = j * size; k < j * size + size; k++) - { - int ag0, ag1; - float agr; - agr = (agp[k] / 360.0f) * (sbin * 2); - ag0 = (int)agr; - ag1 = (ag0 + 1 < sbin * 2) ? ag0 + 1 : 0; - agr = agr - ag0; - cnp[ag0] += (1.0f - agr) * mgp[k] / 255.0f; - cnp[ag1] += agr * mgp[k] / 255.0f; - } - cnp += 2 * sbin; - } - 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.fl; - 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.fl; - cap = ca->data.db; - ccv_zero(db); - float* dbp = db->data.fl; - // 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) \ - { \ - float norm = 1.0f / sqrt(cap[a] + cap[b] + cap[c] + cap[d] + 1e-4f); \ - for (k = 0; k < sbin * 2; k++) \ - { \ - float v = 0.5f * ccv_min(cnp[k] * norm, 0.2f); \ - dbp[5 + sbin + k] += v; \ - dbp[1 + idx] += v; \ - } \ - dbp[sbin * 3 + idx] *= 0.2357f; \ - for (k = 0; k < sbin; k++) \ - { \ - float v = 0.5f * ccv_min((cnp[k] + cnp[k + sbin]) * norm, 0.2f); \ - 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); -} ccv_array_t* ccv_dpm_detect_objects(ccv_dense_matrix_t* a, ccv_dpm_root_classifier_t** _classifier, int count, ccv_dpm_param_t params) { @@ -184,33 +23,12 @@ ccv_array_t* ccv_dpm_detect_objects(ccv_dense_matrix_t* a, ccv_dpm_root_classifi ccv_resample(pyr[0], &pyr[i], 0, (int)(pyr[0]->rows / pow(scale, i)), (int)(pyr[0]->cols / pow(scale, i)), CCV_INTER_AREA); for (i = next; i < scale_upto + next * 2; i++) ccv_sample_down(pyr[i - next], &pyr[i], 0, 0, 0); - for (i = 0; i < 1; i++) // scale_upto + next * 2; i++) + for (i = 0; i < scale_upto + next * 2; i++) { ccv_dense_matrix_t* hog = 0; _ccv_hog(pyr[i], &hog, 0, 9, 8); - _ccv_run_lsvm(hog, ); ccv_matrix_free(hog); } - /* - ccv_dense_matrix_t* hog = 0; - _ccv_hog(pyr[0], &hog, 0, 9, 8); - ccv_dense_matrix_t* b = ccv_dense_matrix_new(pyr[0]->rows, pyr[0]->cols, CCV_8U | CCV_C1, 0, 0); - unsigned char* bptr = b->data.ptr; - for (i = 0; i < b->rows; i++) - { - if (i >= hog->rows * 8) - break; - for (j = 0; j < b->cols; j++) - { - int k = (i / 8) * hog->cols + (j / 8) * 31 + 2; - bptr[j] = ccv_clamp(100 * hog->data.fl[k], 0, 255); - } - bptr += b->step; - } - ccv_serialize(b, "hog.png", 0, CCV_SERIAL_PNG_FILE, 0); - ccv_matrix_free(hog); - ccv_matrix_free(b); - */ if (params.size.height != _classifier[0]->root.size.height * 8 || params.size.width != _classifier[0]->root.size.width * 8) ccv_matrix_free(pyr[0]); for (i = 1; i < scale_upto + next * 2; i++) diff --git a/lib/ccv_numeric.c b/lib/ccv_numeric.c index 8ff8095d5..0e56261f9 100644 --- a/lib/ccv_numeric.c +++ b/lib/ccv_numeric.c @@ -873,3 +873,108 @@ void ccv_convolve_kernel(ccv_dense_matrix_t* x, ccv_convolve_kernel_f func, void #undef for_block ccv_matrix_generate_signature((char*) x->data.ptr, x->rows * x->step, x->sig, 0); } + +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) +{ + assert(!(flag & CCV_L2_NORM) && (flag & CCV_GSEDT)); + ccv_declare_matrix_signature(sig, a->sig != 0, ccv_sign_with_format(64, "ccv_distance_transform(%lf,%lf,%lf,%lf,%d)", dx, dy, dxx, dyy, flag), a->sig, 0); + type = (CCV_GET_DATA_TYPE(type) == CCV_32F) ? CCV_GET_CHANNEL(a->type) | CCV_32F : CCV_GET_CHANNEL(a->type) | CCV_64F; + 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_matrix_return_if_cached(, db); + int i, j, k; + unsigned char* a_ptr = a->data.ptr; + unsigned char* b_ptr = db->data.ptr; + int* v = (int*)alloca(sizeof(int) * ccv_max(db->rows, db->cols)); + double* z = (double*)alloca(sizeof(double) * (ccv_max(db->rows, db->cols) + 1)); + unsigned char* c_ptr = (unsigned char*)alloca(CCV_GET_DATA_TYPE_SIZE(db->type) * db->rows); +#define for_block(_for_set_b, _for_get_b, _for_get_a) \ + if (dxx > 1e-10) \ + { \ + for (i = 0; i < a->rows; i++) \ + { \ + k = 0; \ + v[0] = 0; \ + z[0] = -DBL_MAX; \ + z[1] = DBL_MAX; \ + for (j = 1; j < a->cols; j++) \ + { \ + double s; \ + for (;;) \ + { \ + s = ((_for_get_a(a_ptr, j, 0) + dxx * j * j - dx * j) - (_for_get_a(a_ptr, v[k], 0) + dxx * v[k] * v[k] - dx * v[k])) / (2.0 * dxx * (j - v[k])); \ + if (s > z[k]) break; \ + --k; \ + } \ + ++k; \ + v[k] = j; \ + z[k] = s; \ + z[k + 1] = db->cols * 1e5; \ + } \ + k = 0; \ + for (j = 0; j < a->cols; j++) \ + { \ + while (z[k + 1] < j) \ + ++k; \ + _for_set_b(b_ptr, j, dx * (j - v[k]) + dxx * (j - v[k]) * (j - v[k]) + _for_get_a(a_ptr, v[k], 0), 0); \ + } \ + a_ptr += a->step; \ + b_ptr += db->step; \ + } \ + } else { /* above algorithm cannot handle dxx == 0 properly, below is special casing for that */ \ + for (i = 0; i < a->rows; i++) \ + { \ + for (j = 0; j < a->cols; j++) \ + _for_set_b(b_ptr, j, _for_get_a(a_ptr, j, 0), 0); \ + for (j = 1; j < a->cols; j++) \ + _for_set_b(b_ptr, j, ccv_min(_for_get_b(b_ptr, j, 0), _for_get_b(b_ptr, j - 1, 0) + dx), 0); \ + for (j = a->cols - 2; j >= 0; j--) \ + _for_set_b(b_ptr, j, ccv_min(_for_get_b(b_ptr, j, 0), _for_get_b(b_ptr, j + 1, 0) - dx), 0); \ + a_ptr += a->step; \ + b_ptr += db->step; \ + } \ + } \ + b_ptr = db->data.ptr; \ + if (dyy > 1e-10) \ + { \ + for (j = 0; j < db->cols; j++) \ + { \ + for (i = 0; i < db->rows; i++) \ + _for_set_b(c_ptr, i, _for_get_b(b_ptr + i * db->step, j, 0), 0); \ + k = 0; \ + v[0] = 0; \ + z[0] = -DBL_MAX; \ + z[1] = DBL_MAX; \ + for (i = 1; i < db->rows; i++) \ + { \ + double s; \ + for (;;) \ + { \ + s = ((_for_get_b(c_ptr, i, 0) + dyy * i * i - dy * i) - (_for_get_b(c_ptr, v[k], 0) + dyy * v[k] * v[k] - dy * v[k])) / (2.0 * dyy * (i - v[k])); \ + if (s > z[k]) break; \ + --k; \ + } \ + ++k; \ + v[k] = i; \ + z[k] = s; \ + z[k + 1] = db->rows * 1e5; \ + } \ + k = 0; \ + for (i = 0; i < db->rows; i++) \ + { \ + while (z[k + 1] < i) \ + ++k; \ + _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], 0), 0); \ + } \ + } \ + } else { \ + for (j = 0; j < db->cols; j++) \ + { \ + for (i = 1; i < db->rows; i++) \ + _for_set_b(b_ptr + i * db->step, j, ccv_min(_for_get_b(b_ptr + i * db->step, j, 0), _for_get_b(b_ptr + (i - 1) * db->step, j, 0) + dy), 0); \ + for (i = db->rows - 2; i >= 0; i--) \ + _for_set_b(b_ptr + i * db->step, j, ccv_min(_for_get_b(b_ptr + i * db->step, j, 0), _for_get_b(b_ptr + (i + 1) * db->step, j, 0) - dy), 0); \ + } \ + } + ccv_matrix_setter_getter(db->type, ccv_matrix_getter, a->type, for_block); +#undef for_block +} diff --git a/lib/makefile b/lib/makefile index b8173d6e6..6279734b0 100644 --- a/lib/makefile +++ b/lib/makefile @@ -6,7 +6,7 @@ all: libccv.a clean: rm -f *.o 3rdparty/sha1/*.o 3rdparty/kissfft/*.o libccv.a -libccv.a: ccv_cache.o ccv_memory.o 3rdparty/sha1/sha1.o 3rdparty/kissfft/kiss_fft.o 3rdparty/kissfft/kiss_fftnd.o 3rdparty/kissfft/kiss_fftr.o 3rdparty/kissfft/kiss_fftndr.o ccv_io.o ccv_numeric.o ccv_algebra.o ccv_util.o ccv_basic.o ccv_daisy.o ccv_sift.o ccv_bbf.o ccv_swt.o +libccv.a: ccv_cache.o ccv_memory.o 3rdparty/sha1/sha1.o 3rdparty/kissfft/kiss_fft.o 3rdparty/kissfft/kiss_fftnd.o 3rdparty/kissfft/kiss_fftr.o 3rdparty/kissfft/kiss_fftndr.o ccv_io.o ccv_numeric.o ccv_algebra.o ccv_util.o ccv_basic.o ccv_daisy.o ccv_sift.o ccv_bbf.o ccv_swt.o ccv_dpm.o ar rcs $@ $^ ccv_io.o: ccv_io.c ccv.h io/*.c diff --git a/samples/font.png b/samples/font.png new file mode 100644 index 000000000..fe59399c0 Binary files /dev/null and b/samples/font.png differ diff --git a/samples/geometry.png b/samples/geometry.png new file mode 100644 index 000000000..aebd91003 Binary files /dev/null and b/samples/geometry.png differ diff --git a/samples/pedestrian.png b/samples/pedestrian.png new file mode 100644 index 000000000..00622ca8d Binary files /dev/null and b/samples/pedestrian.png differ diff --git a/samples/street.png b/samples/street.png new file mode 100644 index 000000000..14f49136b Binary files /dev/null and b/samples/street.png differ diff --git a/test/algebra.tests.c b/test/algebra.tests.c index 6739633ea..d30880350 100644 --- a/test/algebra.tests.c +++ b/test/algebra.tests.c @@ -68,7 +68,7 @@ TEST_CASE("summed area table without padding") ptr += dmt->step; } ccv_dense_matrix_t* b = 0; - ccv_sat(dmt, &b, 0, CCV_SAT_NO_PADDING); + ccv_sat(dmt, &b, 0, CCV_NO_PADDING); int sat[60] = { 1, 2, 3, 2, 4, 6, 3, 6, 9, 4, 8, 12, 2, 4, 6, 4, 8, 12, 6, 12, 18, 8, 16, 24, 3, 6, 9, 6, 12, 18, 9, 18, 27, 12, 24, 36, @@ -95,7 +95,7 @@ TEST_CASE("summed area table with padding") ptr += dmt->step; } ccv_dense_matrix_t* b = 0; - ccv_sat(dmt, &b, 0, CCV_SAT_PADDING); + ccv_sat(dmt, &b, 0, CCV_PADDING_ZERO); int sat[72] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 2, 4, 6, 3, 6, 9, 0, 0, 0, 2, 4, 6, 4, 8, 12, 6, 12, 18, diff --git a/test/case_main.h b/test/case_main.h index 598c3f13a..87e65ba4c 100644 --- a/test/case_main.h +++ b/test/case_main.h @@ -259,6 +259,7 @@ int main(int argc, char** argv) if (test_suite->sig == 0x883253372849284B) { printf("\033[0;34m[%d/%d]\033[0;0m \033[1;33m[RUN]\033[0;0m %s ...", j, total, test_suite->name); + fflush(stdout); int result = 0; test_suite->driver(test_suite->name, &result); if (result == 0) diff --git a/test/data/nature.filter.bin b/test/data/nature.filter.bin deleted file mode 100644 index 212a84eb8..000000000 Binary files a/test/data/nature.filter.bin and /dev/null differ diff --git a/test/numeric.tests.c b/test/numeric.tests.c index f810a9aa5..c19b3b7f2 100644 --- a/test/numeric.tests.c +++ b/test/numeric.tests.c @@ -2,6 +2,12 @@ #include "case.h" #include "ccv_case.h" +/* numeric tests are more like functional tests rather than unit tests: + * the following tests contain: + * 1. minimization of the famous rosenbrock function; + * 2. compute ssd with ccv_convolve, and compare the result with naive method + * 3. compare the result from ccv_distance_transform (linear time) with reference implementation from voc-release4 (O(nlog(n))) */ + int rosenbrock(const ccv_dense_matrix_t* x, double* f, ccv_dense_matrix_t* df, void* data) { int* steps = (int*)data; @@ -40,27 +46,175 @@ TEST_CASE("minimize rosenbrock") ccv_matrix_free(x); } -double gaussian(double x, double y, void* data) +static void naive_ssd(ccv_dense_matrix_t* image, ccv_dense_matrix_t* template, ccv_dense_matrix_t* out) { - return exp(-(x * x + y * y) / 20) / sqrt(CCV_PI * 20); + int thw = template->cols / 2; + int thh = template->rows / 2; + int i, j, k, x, y, ch = CCV_GET_CHANNEL(image->type); + unsigned char* i_ptr = image->data.ptr + thh * image->step; + double* o = out->data.db + out->cols * thh; + ccv_zero(out); + for (i = thh; i < image->rows - thh - 1; i++) + { + for (j = thw; j < image->cols - thw - 1; j++) + { + unsigned char* t_ptr = template->data.ptr; + unsigned char* j_ptr = i_ptr - thh * image->step; + o[j] = 0; + for (y = -thh; y <= thh; y++) + { + for (x = -thw; x <= thw; x++) + for (k = 0; k < ch; k++) + o[j] += (j_ptr[(x + j) * ch + k] - t_ptr[(x + thw) * ch + k]) * (j_ptr[(x + j) * ch + k] - t_ptr[(x + thw) * ch + k]); + t_ptr += template->step; + j_ptr += image->step; + } + } + i_ptr += image->step; + o += out->cols; + } } -TEST_CASE("convolution on Gaussian kernel") +TEST_CASE("convolution ssd (sum of squared differences) v.s. naive ssd") { - ccv_dense_matrix_t* image = 0; - ccv_read("../samples/nature.png", &image, CCV_IO_GRAY | CCV_IO_ANY_FILE); - ccv_dense_matrix_t* gray = 0;ccv_dense_matrix_new(image->rows, image->cols, CCV_32F | CCV_C1, 0, 0); - ccv_shift(image, (ccv_matrix_t**)&gray, CCV_32F | CCV_C1, 0, 0); - ccv_dense_matrix_t* kernel = ccv_dense_matrix_new(10, 10, CCV_32F | CCV_C1, 0, 0); - ccv_filter_kernel(kernel, gaussian, 0); - ccv_normalize(kernel, (ccv_matrix_t**)&kernel, 0, CCV_L1_NORM); - ccv_dense_matrix_t* x = 0; - ccv_filter(gray, kernel, (ccv_matrix_t**)&x, 0); - REQUIRE_MATRIX_FILE_EQ(x, "data/nature.filter.bin", "should apply Gaussian filter with FFTW on nature.png with sigma = sqrt(10)"); - ccv_matrix_free(image); - ccv_matrix_free(gray); - ccv_matrix_free(kernel); - ccv_matrix_free(x); + ccv_dense_matrix_t* street = 0; + ccv_dense_matrix_t* pedestrian = 0; + ccv_read("../samples/pedestrian.png", &pedestrian, CCV_IO_ANY_FILE); + ccv_read("../samples/street.png", &street, CCV_IO_ANY_FILE); + ccv_dense_matrix_t* result = 0; + ccv_convolve(street, pedestrian, &result, CCV_64F, 0); + ccv_dense_matrix_t* square = 0; + ccv_multiply(street, street, (ccv_matrix_t**)&square, 0); + ccv_dense_matrix_t* sat = 0; + ccv_sat(square, &sat, 0, CCV_PADDING_ZERO); + ccv_matrix_free(square); + double sum[] = {0, 0, 0}; + int i, j, k; + int ch = CCV_GET_CHANNEL(street->type); + unsigned char* p_ptr = pedestrian->data.ptr; +#define for_block(_, _for_get) \ + for (i = 0; i < pedestrian->rows; i++) \ + { \ + for (j = 0; j < pedestrian->cols; j++) \ + for (k = 0; k < ch; k++) \ + sum[k] += _for_get(p_ptr, j * ch + k, 0) * _for_get(p_ptr, j * ch + k, 0); \ + p_ptr += pedestrian->step; \ + } + ccv_matrix_getter(pedestrian->type, for_block); +#undef for_block + int phw = pedestrian->cols / 2; + int phh = pedestrian->rows / 2; + ccv_dense_matrix_t* output = ccv_dense_matrix_new(street->rows, street->cols, CCV_64F | CCV_C1, 0, 0); + ccv_zero(output); + unsigned char* s_ptr = sat->data.ptr + sat->step * phh; + unsigned char* r_ptr = result->data.ptr + result->step * phh; + double* o_ptr = output->data.db + output->cols * phh; +#define for_block(_for_get_s, _for_get_r) \ + for (i = phh; i < output->rows - phh - 1; i++) \ + { \ + for (j = phw; j < output->cols - phw - 1; j++) \ + { \ + o_ptr[j] = 0; \ + for (k = 0; k < ch; k++) \ + { \ + o_ptr[j] += (_for_get_s(s_ptr + sat->step * ccv_min(phh + 1, sat->rows - i - 1), ccv_min(j + phw + 1, sat->cols - 1) * ch + k, 0) \ + - _for_get_s(s_ptr + sat->step * ccv_min(phh + 1, sat->rows - i - 1), ccv_max(j - phw, 0) * ch + k, 0) \ + + _for_get_s(s_ptr + sat->step * ccv_max(-phh, -i), ccv_max(j - phw, 0) * ch + k, 0) \ + - _for_get_s(s_ptr + sat->step * ccv_max(-phh, -i), ccv_min(j + phw + 1, sat->cols - 1) * ch + k, 0)) \ + + sum[k] - 2.0 * _for_get_r(r_ptr, j * ch + k, 0); \ + } \ + } \ + s_ptr += sat->step; \ + r_ptr += result->step; \ + o_ptr += output->cols; \ + } + ccv_matrix_getter(sat->type, ccv_matrix_getter_a, result->type, for_block); +#undef for_block + ccv_matrix_free(result); + ccv_matrix_free(sat); + ccv_dense_matrix_t* final = 0; + ccv_slice(output, (ccv_matrix_t**)&final, 0, phh, phw, output->rows - phh * 2, output->cols - phw * 2); + ccv_zero(output); + naive_ssd(street, pedestrian, output); + ccv_dense_matrix_t* ref = 0; + ccv_slice(output, (ccv_matrix_t**)&ref, 0, phh, phw, output->rows - phh * 2, output->cols - phw * 2); + ccv_matrix_free(output); + ccv_matrix_free(pedestrian); + ccv_matrix_free(street); + REQUIRE_MATRIX_EQ(ref, final, "ssd computed by convolution doesn't match the one computed by naive method"); + ccv_matrix_free(final); + ccv_matrix_free(ref); +} + +// divide & conquer method for distance transform (copied directly from dpm-matlab (voc-release4) + +static inline int square(int x) { return x*x; } + +// dt helper function +void dt_helper(double *src, double *dst, int *ptr, int step, + int s1, int s2, int d1, int d2, double a, double b) { + if (d2 >= d1) { + int d = (d1+d2) >> 1; + int s = s1; + for (int p = s1+1; p <= s2; p++) + if (src[s*step] + a*square(d-s) + b*(d-s) > + src[p*step] + a*square(d-p) + b*(d-p)) + s = p; + dst[d*step] = src[s*step] + a*square(d-s) + b*(d-s); + ptr[d*step] = s; + dt_helper(src, dst, ptr, step, s1, s, d1, d-1, a, b); + dt_helper(src, dst, ptr, step, s, s2, d+1, d2, a, b); + } +} + +// dt of 1d array +void dt1d(double *src, double *dst, int *ptr, int step, int n, + double a, double b) { + dt_helper(src, dst, ptr, step, 0, n-1, 0, n-1, a, b); +} + +void daq_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, double dx, double dy, double dxx, double dyy) +{ + ccv_dense_matrix_t* dc = ccv_dense_matrix_new(a->rows, a->cols, CCV_64F | CCV_C1, 0, 0); + ccv_dense_matrix_t* db = *b = ccv_dense_matrix_new(a->rows, a->cols, CCV_64F | CCV_C1, 0, 0); + unsigned char* a_ptr = a->data.ptr; + double* b_ptr = db->data.db; + int i, j; +#define for_block(_, _for_get) \ + for (i = 0; i < a->rows; i++) \ + { \ + for (j = 0; j < a->cols; j++) \ + b_ptr[j] = _for_get(a_ptr, j, 0); \ + b_ptr += db->cols; \ + a_ptr += a->step; \ + } + ccv_matrix_getter(a->type, for_block); +#undef for_block + int* ix = (int*)calloc(a->cols * a->rows, sizeof(int)); + int* iy = (int*)calloc(a->cols * a->rows, sizeof(int)); + b_ptr = db->data.db; + double* c_ptr = dc->data.db; + for (i = 0; i < a->rows; i++) + dt1d(b_ptr + i * a->cols, c_ptr + i * a->cols, ix + i * a->cols, 1, a->cols, dxx, dx); + for (j = 0; j < a->cols; j++) + dt1d(c_ptr + j, b_ptr + j, iy + j, a->cols, a->rows, dyy, dy); + free(ix); + free(iy); + ccv_matrix_free(dc); +} + +TEST_CASE("ccv_distance_transform (linear time) v.s. distance transform using divide & conquer (O(nlog(n)))") +{ + ccv_dense_matrix_t* geometry = 0; + ccv_read("../samples/geometry.png", &geometry, CCV_IO_GRAY | CCV_IO_ANY_FILE); + ccv_dense_matrix_t* distance = 0; + ccv_distance_transform(geometry, &distance, 0, 1, 1, 0.1, 0.1, CCV_GSEDT); + ccv_dense_matrix_t* ref = 0; + daq_distance_transform(geometry, &ref, 1, 1, 0.1, 0.1); + ccv_matrix_free(geometry); + REQUIRE_MATRIX_EQ(distance, ref, "distance transform computed by ccv_distance_transform doesn't match the one computed by divide & conquer (voc-release4)"); + ccv_matrix_free(ref); + ccv_matrix_free(distance); } #include "case_main.h"