fixed unit test which single-precision broke

pteichman · Apr 5, 2012 · f6bc465 · f6bc465
1 parent d2fb75e
commit f6bc465
Show file tree

Hide file tree

Showing 6 changed files with 65 additions and 54 deletions.
diff --git a/lib/ccv_basic.c b/lib/ccv_basic.c
@@ -1067,7 +1067,7 @@ void ccv_blur(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, double si
 	ccv_matrix_return_if_cached(, db);
 	int fsz = ccv_max(1, (int)(4.0 * sigma + 1.0 - 1e-8)) * 2 + 1;
 	int hfz = fsz / 2;
-	unsigned char* buf = (unsigned char*)alloca(sizeof(double) * (fsz + ccv_max(a->rows, a->cols * CCV_GET_CHANNEL(a->type))));
+	unsigned char* buf = (unsigned char*)alloca(sizeof(double) * ccv_max(fsz + a->rows, (fsz + 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);

diff --git a/lib/ccv_numeric.c b/lib/ccv_numeric.c
@@ -476,7 +476,7 @@ static int _ccv_get_optimal_fft_size(int size)
 static void _ccv_filter_fftw(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t* d, int padding_pattern)
 {
 	int ch = CCV_GET_CHANNEL(a->type);
-	int fft_type = (CCV_GET_DATA_TYPE(b->type) == CCV_8U || CCV_GET_DATA_TYPE(b->type) == CCV_32F) ? CCV_32F : CCV_64F;
+	int fft_type = (CCV_GET_DATA_TYPE(d->type) == CCV_8U || CCV_GET_DATA_TYPE(d->type) == CCV_32F) ? CCV_32F : CCV_64F;
 	int rows = ccv_min(a->rows + b->rows - 1, _ccv_get_optimal_fft_size(b->rows * 3));
 	int cols = ccv_min(a->cols + b->cols - 1, _ccv_get_optimal_fft_size(b->cols * 3));
 	int cols_2c = 2 * (cols / 2 + 1);
@@ -654,7 +654,7 @@ static void _ccv_filter_fftw(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_d
 static void _ccv_filter_kissfft(ccv_dense_matrix_t* a, ccv_dense_matrix_t* b, ccv_dense_matrix_t* d, int padding_pattern)
 {
 	int ch = CCV_GET_CHANNEL(a->type);
-	int fft_type = (CCV_GET_DATA_TYPE(b->type) == CCV_8U || CCV_GET_DATA_TYPE(b->type) == CCV_32F) ? CCV_32F : CCV_64F;
+	int fft_type = (CCV_GET_DATA_TYPE(d->type) == CCV_8U || CCV_GET_DATA_TYPE(d->type) == CCV_32F) ? CCV_32F : CCV_64F;
 	int rows = ((ccv_min(a->rows + b->rows - 1, kiss_fftr_next_fast_size_real(b->rows * 3)) + 1) >> 1) << 1;
 	int cols = ((ccv_min(a->cols + b->cols - 1, kiss_fftr_next_fast_size_real(b->cols * 3)) + 1) >> 1) << 1;
 	int ndim[] = {rows, cols};
@@ -988,25 +988,26 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 	unsigned char* a_ptr = a->data.u8;
 	unsigned char* b_ptr = db->data.u8;
 	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);
 	int* x_ptr = mx ? mx->data.i32 : 0;
 	int* y_ptr = my ? my->data.i32 : 0;
-#define for_block(_for_set_b, _for_get_b, _for_get_a) \
-	if (dxx > 1e-10) \
+#define for_block(_for_max, _for_type_b, _for_set_b, _for_get_b, _for_get_a) \
+	_for_type_b _dx = dx, _dy = dy, _dxx = dxx, _dyy = dyy; \
+	_for_type_b* z = (_for_type_b*)alloca(sizeof(_for_type_b) * (ccv_max(db->rows, db->cols) + 1)); \
+	if (_dxx > 1e-6) \
 	{ \
 		for (i = 0; i < a->rows; i++) \
 		{ \
 			k = 0; \
 			v[0] = 0; \
-			z[0] = -DBL_MAX; \
-			z[1] = DBL_MAX; \
+			z[0] = -_for_max; \
+			z[1] = _for_max; \
 			for (j = 1; j < a->cols; j++) \
 			{ \
-				double s; \
+				_for_type_b s; \
 				for (;;) \
 				{ \
-					s = ((SGN _for_get_a(a_ptr, j, 0) + dxx * j * j - dx * j) - (SGN _for_get_a(a_ptr, v[k], 0) + dxx * v[k] * v[k] - dx * v[k])) / (2.0 * dxx * (j - v[k])); \
+					s = ((SGN _for_get_a(a_ptr, j, 0) + _dxx * j * j - _dx * j) - (SGN _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; \
 				} \
@@ -1022,7 +1023,7 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 				{ \
 					while (z[k + 1] < j) \
 						++k; \
-					_for_set_b(b_ptr, j, dx * (j - v[k]) + dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k], 0), 0); \
+					_for_set_b(b_ptr, j, _dx * (j - v[k]) + _dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k], 0), 0); \
 					x_ptr[j] = j - v[k]; \
 				} \
 				x_ptr += mx->cols; \
@@ -1031,7 +1032,7 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 				{ \
 					while (z[k + 1] < j) \
 						++k; \
-					_for_set_b(b_ptr, j, dx * (j - v[k]) + dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k], 0), 0); \
+					_for_set_b(b_ptr, j, _dx * (j - v[k]) + _dxx * (j - v[k]) * (j - v[k]) SGN _for_get_a(a_ptr, v[k], 0), 0); \
 				} \
 			} \
 			a_ptr += a->step; \
@@ -1044,30 +1045,30 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 			for (j = 0; j < a->cols; j++) \
 				_for_set_b(b_ptr, j, SGN _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_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); \
+				_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.u8; \
-	if (dyy > 1e-10) \
+	if (_dyy > 1e-6) \
 	{ \
 		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; \
+			z[0] = -_for_max; \
+			z[1] = _for_max; \
 			for (i = 1; i < db->rows; i++) \
 			{ \
-				double s; \
+				_for_type_b 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])); \
+					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; \
 				} \
@@ -1083,7 +1084,7 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 				{ \
 					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); \
+					_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); \
 					y_ptr[i * my->cols] = i - v[k]; \
 				} \
 				++y_ptr; \
@@ -1092,7 +1093,7 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 				{ \
 					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); \
+					_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); \
 				} \
 			} \
 		} \
@@ -1101,19 +1102,29 @@ void ccv_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int t
 		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_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); \
+				_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); \
 		} \
 	}
 	if (flag & CCV_NEGATIVE)
 	{
 #define SGN -
-		ccv_matrix_setter_getter(db->type, ccv_matrix_getter, a->type, for_block);
+		if (db->type & CCV_64F)
+		{
+			ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, DBL_MAX);
+		} else {
+			ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, FLT_MAX);
+		}
 #undef SGN
 	} else {
 #define SGN +
-		ccv_matrix_setter_getter(db->type, ccv_matrix_getter, a->type, for_block);
+		if (db->type & CCV_64F)
+		{
+			ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, DBL_MAX);
+		} else {
+			ccv_matrix_typeof_setter_getter(db->type, ccv_matrix_getter, a->type, for_block, FLT_MAX);
+		}
 #undef SGN
 	}
 #undef for_block

diff --git a/lib/ccv_util.c b/lib/ccv_util.c
@@ -457,7 +457,7 @@ int ccv_matrix_eq(ccv_matrix_t* a, ccv_matrix_t* b)
 		{ \
 			for (j = 0; j < da->cols * ch; j++) \
 			{ \
-				if (fabs(_for_get(b_ptr, j, 0) - _for_get(a_ptr, j, 0)) > 1e-6) \
+				if (fabs(_for_get(b_ptr, j, 0) - _for_get(a_ptr, j, 0)) > 1e-4) \
 					return -1; \
 			} \
 			a_ptr += da->step; \

diff --git a/lib/makefile b/lib/makefile
@@ -1,5 +1,5 @@
 CC = clang
-CFLAGS = -O3 -msse2 -Wall -D HAVE_LIBJPEG -D HAVE_CBLAS -D HAVE_GSL -D HAVE_LIBPNG -D HAVE_FFTW3 # -fopenmp -D USE_OPENMP
+CFLAGS = -O3 -msse2 -ffast-math -mtune=native -Wall -D HAVE_LIBJPEG -D HAVE_CBLAS -D HAVE_GSL -D HAVE_LIBPNG -D HAVE_FFTW3 # -fopenmp -D USE_OPENMP
 
 all: libccv.a
 

diff --git a/test/data/nature.blur.bin b/test/data/nature.blur.bin
diff --git a/test/numeric.tests.c b/test/numeric.tests.c
@@ -181,8 +181,8 @@ TEST_CASE("convolution ssd (sum of squared differences) v.s. naive ssd")
 static inline int square(int x) { return x*x; }
 
 // dt helper function
-void dt_min_helper(double *src, double *dst, int *ptr, int step, 
-	       int s1, int s2, int d1, int d2, double a, double b) {
+void dt_min_helper(float *src, float *dst, int *ptr, int step, 
+	       int s1, int s2, int d1, int d2, float a, float b) {
  if (d2 >= d1) {
    int d = (d1+d2) >> 1;
    int s = s1;
@@ -198,17 +198,17 @@ void dt_min_helper(double *src, double *dst, int *ptr, int step,
 }
 
 // dt of 1d array
-void dt_min1d(double *src, double *dst, int *ptr, int step, int n, 
-	  double a, double b) {
+void dt_min1d(float *src, float *dst, int *ptr, int step, int n, 
+	  float a, float b) {
   dt_min_helper(src, dst, ptr, step, 0, n-1, 0, n-1, a, b);
 }
 
 void daq_min_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);
+	ccv_dense_matrix_t* dc = ccv_dense_matrix_new(a->rows, a->cols, CCV_32F | CCV_C1, 0, 0);
+	ccv_dense_matrix_t* db = *b = ccv_dense_matrix_new(a->rows, a->cols, CCV_32F | CCV_C1, 0, 0);
 	unsigned char* a_ptr = a->data.u8;
-	double* b_ptr = db->data.f64;
+	float* b_ptr = db->data.f32;
 	int i, j;
 #define for_block(_, _for_get) \
 	for (i = 0; i < a->rows; i++) \
@@ -222,8 +222,8 @@ void daq_min_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, d
 #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.f64;
-	double* c_ptr = dc->data.f64;
+	b_ptr = db->data.f32;
+	float* c_ptr = dc->data.f32;
 	for (i = 0; i < a->rows; i++)
 		dt_min1d(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++)
@@ -238,10 +238,10 @@ TEST_CASE("ccv_distance_transform (linear time) v.s. distance transform using di
 	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;
-	double dx = 1;
-	double dy = 1;
-	double dxx = 0.4;
-	double dyy = 0.4;
+	double dx = 0;
+	double dy = 0;
+	double dxx = 1;
+	double dyy = 1;
 	ccv_distance_transform(geometry, &distance, 0, 0, 0, 0, 0, dx, dy, dxx, dyy, CCV_GSEDT);
 	ccv_dense_matrix_t* ref = 0;
 	daq_min_distance_transform(geometry, &ref, dx, dy, dxx, dyy);
@@ -252,8 +252,8 @@ TEST_CASE("ccv_distance_transform (linear time) v.s. distance transform using di
 }
 
 // dt helper function
-void dt_max_helper(double *src, double *dst, int *ptr, int step, 
-	       int s1, int s2, int d1, int d2, double a, double b) {
+void dt_max_helper(float *src, float *dst, int *ptr, int step, 
+	       int s1, int s2, int d1, int d2, float a, float b) {
  if (d2 >= d1) {
    int d = (d1+d2) >> 1;
    int s = s1;
@@ -269,17 +269,17 @@ void dt_max_helper(double *src, double *dst, int *ptr, int step,
 }
 
 // dt of 1d array
-void dt_max1d(double *src, double *dst, int *ptr, int step, int n, 
-	  double a, double b) {
+void dt_max1d(float *src, float *dst, int *ptr, int step, int n, 
+	  float a, float b) {
   dt_max_helper(src, dst, ptr, step, 0, n-1, 0, n-1, a, b);
 }
 
 void daq_max_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);
+	ccv_dense_matrix_t* dc = ccv_dense_matrix_new(a->rows, a->cols, CCV_32F | CCV_C1, 0, 0);
+	ccv_dense_matrix_t* db = *b = ccv_dense_matrix_new(a->rows, a->cols, CCV_32F | CCV_C1, 0, 0);
 	unsigned char* a_ptr = a->data.u8;
-	double* b_ptr = db->data.f64;
+	float* b_ptr = db->data.f32;
 	int i, j;
 #define for_block(_, _for_get) \
 	for (i = 0; i < a->rows; i++) \
@@ -293,8 +293,8 @@ void daq_max_distance_transform(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, d
 #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.f64;
-	double* c_ptr = dc->data.f64;
+	b_ptr = db->data.f32;
+	float* c_ptr = dc->data.f32;
 	for (i = 0; i < a->rows; i++)
 		dt_max1d(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++)
@@ -309,17 +309,17 @@ TEST_CASE("ccv_distance_transform to compute max distance")
 	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;
-	double dx = 0;
-	double dy = 0;
-	double dxx = 1;
-	double dyy = 1;
+	double dx = 1;
+	double dy = 1;
+	double dxx = 0.4;
+	double dyy = 0.4;
 	ccv_distance_transform(geometry, &distance, 0, 0, 0, 0, 0, dx, dy, dxx, dyy, CCV_NEGATIVE | CCV_GSEDT);
 	ccv_dense_matrix_t* ref = 0;
 	daq_max_distance_transform(geometry, &ref, dx, dy, dxx, dyy);
 	ccv_matrix_free(geometry);
 	int i;
 	for (i = 0; i < distance->rows * distance->cols; i++)
-		distance->data.f64[i] = -distance->data.f64[i];
+		distance->data.f32[i] = -distance->data.f32[i];
 	REQUIRE_MATRIX_EQ(distance, ref, "maximum distance transform computed by negate ccv_distance_transform doesn't match the one computed by divide & conquer");
 	ccv_matrix_free(ref);
 	ccv_matrix_free(distance);