Merge pull request #55 from claudiofelber/master

Replaced existing threshold function with Otsu threshold algorithm

lib/identify.c

@@ -174,60 +174,55 @@ static void flood_fill_seed(struct quirc *q, int x, int y, int from, int to,
  * Adaptive thresholding
  */
 
-#define THRESHOLD_S_MIN		1
+uint8_t otsu(const struct quirc *q)
-#define THRESHOLD_S_DEN		8
-#define THRESHOLD_T		5
-
-static void threshold(struct quirc *q)
 {
-	int x, y;
+	int numPixels = q->w * q->h;
-	int avg_w = 0;
-	int avg_u = 0;
-	int threshold_s = q->w / THRESHOLD_S_DEN;
-	quirc_pixel_t *row = q->pixels;
 
-	/*
+	// Calculate histogram
-	 * Ensure a sane, non-zero value for threshold_s.
+	const int HISTOGRAM_SIZE = 256;
-	 *
+	unsigned int histogram[HISTOGRAM_SIZE];
-	 * threshold_s can be zero if the image width is small. We need to avoid
+	memset(histogram, 0, (HISTOGRAM_SIZE) * sizeof(unsigned int));
-	 * SIGFPE as it will be used as divisor.
+	uint8_t* ptr = q->image;
-	 */
+	int length = numPixels;
-	if (threshold_s < THRESHOLD_S_MIN)
+	while (length--) {
-		threshold_s = THRESHOLD_S_MIN;
+		uint8_t value = *ptr++;
-
+		histogram[value]++;
-	for (y = 0; y < q->h; y++) {
-		memset(q->row_average, 0, q->w * sizeof(int));
-
-		for (x = 0; x < q->w; x++) {
-			int w, u;
-
-			if (y & 1) {
-				w = x;
-				u = q->w - 1 - x;
-			} else {
-				w = q->w - 1 - x;
-				u = x;
-			}
-
-			avg_w = (avg_w * (threshold_s - 1)) /
-				threshold_s + row[w];
-			avg_u = (avg_u * (threshold_s - 1)) /
-				threshold_s + row[u];
-
-			q->row_average[w] += avg_w;
-			q->row_average[u] += avg_u;
-		}
-
-		for (x = 0; x < q->w; x++) {
-			if (row[x] < q->row_average[x] *
-			    (100 - THRESHOLD_T) / (200 * threshold_s))
-				row[x] = QUIRC_PIXEL_BLACK;
-			else
-				row[x] = QUIRC_PIXEL_WHITE;
-		}
-
-		row += q->w;
 	}
+
+	// Calculate weighted sum of histogram values
+	int sum = 0;
+	for (int i = 0; i < HISTOGRAM_SIZE; ++i) {
+		sum += i * histogram[i];
+	}
+
+	// Compute threshold
+	int sumB = 0;
+	int q1 = 0;
+	double max = 0;
+	uint8_t threshold = 0;
+	for (int i = 0; i < HISTOGRAM_SIZE; ++i) {
+		// Weighted background
+		q1 += histogram[i];
+		if (q1 == 0)
+			continue;
+
+		// Weighted foreground
+		const int q2 = numPixels - q1;
+		if (q2 == 0)
+			break;
+
+		sumB += i * histogram[i];
+		const double m1 = (double)sumB / q1;
+		const double m2 = ((double)sum - sumB) / q2;
+		const double m1m2 = m1 - m2;
+		const double variance = m1m2 * m1m2 * q1 * q2;
+		if (variance > max) {
+			threshold = i;
+			max = variance;
+		}
+	}
+
+	return threshold;
 }
 
 static void area_count(void *user_data, int y, int left, int right)
@@ -1074,17 +1069,18 @@ static void test_grouping(struct quirc *q, int i)
 	test_neighbours(q, i, &hlist, &vlist);
 }
 
-static void pixels_setup(struct quirc *q)
+static void pixels_setup(struct quirc *q, uint8_t threshold)
 {
 	if (sizeof(*q->image) == sizeof(*q->pixels)) {
 		q->pixels = (quirc_pixel_t *)q->image;
-	} else {
+	}
-		int x, y;
+
-		for (y = 0; y < q->h; y++) {
+	uint8_t* source = q->image;
-			for (x = 0; x < q->w; x++) {
+	quirc_pixel_t* dest = q->pixels;
-				q->pixels[y * q->w + x] = q->image[y * q->w + x];
+	int length = q->w * q->h;
-			}
+	while (length--) {
-		}
+		uint8_t value = *source++;
+		*dest++ = (value < threshold) ? QUIRC_PIXEL_BLACK : QUIRC_PIXEL_WHITE;
 	}
 }
 
@@ -1106,8 +1102,8 @@ void quirc_end(struct quirc *q)
 {
 	int i;
 
-	pixels_setup(q);
+	uint8_t threshold = otsu(q);
-	threshold(q);
+	pixels_setup(q, threshold);
 
 	for (i = 0; i < q->h; i++)
 		finder_scan(q, i);

lib/quirc.c

@@ -41,7 +41,6 @@ void quirc_destroy(struct quirc *q)
 	   same size, so we need to be careful here to avoid a double free */
 	if (sizeof(*q->image) != sizeof(*q->pixels))
 		free(q->pixels);
-	free(q->row_average);
 	free(q);
 }
 
@@ -49,7 +48,6 @@ int quirc_resize(struct quirc *q, int w, int h)
 {
 	uint8_t		*image  = NULL;
 	quirc_pixel_t	*pixels = NULL;
-	int		*row_average = NULL;
 
 	/*
 	 * XXX: w and h should be size_t (or at least unsigned) as negatives
@@ -88,11 +86,6 @@ int quirc_resize(struct quirc *q, int w, int h)
 			goto fail;
 	}
 
-	/* alloc a new buffer for q->row_average */
-	row_average = calloc(w, sizeof(int));
-	if (!row_average)
-		goto fail;
-
 	/* alloc succeeded, update `q` with the new size and buffers */
 	q->w = w;
 	q->h = h;
@@ -102,15 +95,12 @@ int quirc_resize(struct quirc *q, int w, int h)
 		free(q->pixels);
 		q->pixels = pixels;
 	}
-	free(q->row_average);
-	q->row_average = row_average;
 
 	return 0;
 	/* NOTREACHED */
 fail:
 	free(image);
 	free(pixels);
-	free(row_average);
 
 	return -1;
 }

lib/quirc_internal.h

@@ -77,7 +77,6 @@ struct quirc_grid {
 struct quirc {
 	uint8_t			*image;
 	quirc_pixel_t		*pixels;
-	int			*row_average; /* used by threshold() */
 	int			w;
 	int			h;