Add an option to use single-precision floating point in calculations

On some platforms, such as microcontrollers with a single-precision
FPU, operations on double type can be significantly slower compared
to the float type. For example, on ESP32-S3 microcontroller, decoding
a QR code of a certain size may take 1700ms when 'double' is used and
just 250ms when 'float' is used.

This commit adds two options to allow for such optimizations:

- QUIRC_FLOAT_TYPE: if defined, it is the type name to use in floating
  point calculations. Can be set, for example, using:
    CFLAGS += -DQUIRC_FLOAT_TYPE=float

- QUIRC_USE_TGMATH: if defined, Quirc will internally use <tgmath.h>
  header, instead of <math.h>. This C99-or-later header allows the
  program to call type-generic functions, such as 'sqrt', and the
  calls will be dispatched to the correct implementation (sqrtf, sqrt,
  sqrtl) depending on the actual argument type. Without setting this
  option, the benefit of -DQUIRC_FLOAT_TYPE=float would be limited as
  the double-precision versions of math functions would still be used.

The change is backwards compatible with existing applications. If
these macros are not defined, the behavior is the same as before.

README.md

@@ -102,10 +102,10 @@ the default macros like CFLAGS from sys.mk can cause unintended effects.
 
 Installation
 ------------
-To build the library and associated demos/tests, type `make`. If you need to
+To build the library and associated demos/tests, type `make`.
-decode "large" image files build with `CFLAGS="-DQUIRC_MAX_REGIONS=65534" make`
+
-instead. Note that this will increase the memory usage, it is discouraged for
+Several options can be adjusted at compile time by passing additional arguments
-low resource devices (i.e. embedded).
+to `make`. See [Compile-time options](#compile-time-options) section below for details.
 
 Type `make install` to install the library, header file and camera demos.
 
@@ -236,6 +236,34 @@ decode attempt with the flipped image data whenever you get an ECC failure:
         printf("Data: %s\n", data.payload);
 ```
 
+Compile-time options
+--------------------
+
+The following compile-time options can be used to adjust the library to a
+particular use case.
+
+Each option is a C preprocessor macro. To set an option, add it to CFLAGS
+using `-DOPTION=VALUE` syntax, for example:
+```bash
+make CFLAGS="-DQUIRC_MAX_REGIONS=65534"
+```
+
+* `QUIRC_MAX_REGIONS`: If you need to decode "large" image files, set
+   `QUIRC_MAX_REGIONS=65534`. Note that since this will increase the memory
+   usage, it is discouraged for low resource devices (i.e. embedded).
+
+* `QUIRC_FLOAT_TYPE`: If defined, it sets the type name to use
+   in floating point calculations. For example, on an embedded system
+   with only a single precision FPU, set `QUIRC_FLOAT_TYPE=float` to
+   improve performance.
+
+* `QUIRC_USE_TGMATH`: if defined, quirc will internally use `<tgmath.h>`
+   header instead of `<math.h>`, ensuring that the math function calls
+   use the same precision as the arguments. Define this option if you are
+   setting `QUIRC_FLOAT_TYPE=float` and the compiler supports C99 or later
+   language standard. 
+
+
 Copyright
 ---------
 Copyright (C) 2010-2012 Daniel Beer <<dlbeer@gmail.com>>

lib/identify.c

@@ -18,7 +18,11 @@
 #include <string.h>
 #include <stdbool.h>
 #include <stdlib.h>
+#ifdef QUIRC_USE_TGMATH
+#include <tgmath.h>
+#else
 #include <math.h>
+#endif // QUIRC_USE_TGMATH
 #include "quirc_internal.h"
 
 /************************************************************************
@@ -62,21 +66,21 @@ static int line_intersect(const struct quirc_point *p0,
 	return 1;
 }
 
-static void perspective_setup(double *c,
+static void perspective_setup(quirc_float_t *c,
 			      const struct quirc_point *rect,
-			      double w, double h)
+			      quirc_float_t w, quirc_float_t h)
 {
-	double x0 = rect[0].x;
+	quirc_float_t x0 = rect[0].x;
-	double y0 = rect[0].y;
+	quirc_float_t y0 = rect[0].y;
-	double x1 = rect[1].x;
+	quirc_float_t x1 = rect[1].x;
-	double y1 = rect[1].y;
+	quirc_float_t y1 = rect[1].y;
-	double x2 = rect[2].x;
+	quirc_float_t x2 = rect[2].x;
-	double y2 = rect[2].y;
+	quirc_float_t y2 = rect[2].y;
-	double x3 = rect[3].x;
+	quirc_float_t x3 = rect[3].x;
-	double y3 = rect[3].y;
+	quirc_float_t y3 = rect[3].y;
 
-	double wden = w * (x2*y3 - x3*y2 + (x3-x2)*y1 + x1*(y2-y3));
+	quirc_float_t wden = w * (x2*y3 - x3*y2 + (x3-x2)*y1 + x1*(y2-y3));
-	double hden = h * (x2*y3 + x1*(y2-y3) - x3*y2 + (x3-x2)*y1);
+	quirc_float_t hden = h * (x2*y3 + x1*(y2-y3) - x3*y2 + (x3-x2)*y1);
 
 	c[0] = (x1*(x2*y3-x3*y2) + x0*(-x2*y3+x3*y2+(x2-x3)*y1) +
 		x1*(x3-x2)*y0) / wden;
@@ -93,24 +97,24 @@ static void perspective_setup(double *c,
 		hden;
 }
 
-static void perspective_map(const double *c,
+static void perspective_map(const quirc_float_t *c,
-			    double u, double v, struct quirc_point *ret)
+			    quirc_float_t u, quirc_float_t v, struct quirc_point *ret)
 {
-	double den = c[6]*u + c[7]*v + 1.0;
+	quirc_float_t den = c[6]*u + c[7]*v + 1.0;
-	double x = (c[0]*u + c[1]*v + c[2]) / den;
+	quirc_float_t x = (c[0]*u + c[1]*v + c[2]) / den;
-	double y = (c[3]*u + c[4]*v + c[5]) / den;
+	quirc_float_t y = (c[3]*u + c[4]*v + c[5]) / den;
 
 	ret->x = (int) rint(x);
 	ret->y = (int) rint(y);
 }
 
-static void perspective_unmap(const double *c,
+static void perspective_unmap(const quirc_float_t *c,
 			      const struct quirc_point *in,
-			      double *u, double *v)
+			      quirc_float_t *u, quirc_float_t *v)
 {
-	double x = in->x;
+	quirc_float_t x = in->x;
-	double y = in->y;
+	quirc_float_t y = in->y;
-	double den = -c[0]*c[7]*y + c[1]*c[6]*y + (c[3]*c[7]-c[4]*c[6])*x +
+	quirc_float_t den = -c[0]*c[7]*y + c[1]*c[6]*y + (c[3]*c[7]-c[4]*c[6])*x +
 		c[0]*c[4] - c[1]*c[3];
 
 	*u = -(c[1]*(y-c[5]) - c[2]*c[7]*y + (c[5]*c[7]-c[4])*x + c[2]*c[4]) /
@@ -299,16 +303,16 @@ static uint8_t otsu(const struct quirc *q)
 	}
 
 	// Calculate weighted sum of histogram values
-	double sum = 0;
+	quirc_float_t sum = 0;
 	unsigned int i = 0;
 	for (i = 0; i <= UINT8_MAX; ++i) {
 		sum += i * histogram[i];
 	}
 
 	// Compute threshold
-	double sumB = 0;
+	quirc_float_t sumB = 0;
 	unsigned int q1 = 0;
-	double max = 0;
+	quirc_float_t max = 0;
 	uint8_t threshold = 0;
 	for (i = 0; i <= UINT8_MAX; ++i) {
 		// Weighted background
@@ -322,10 +326,10 @@ static uint8_t otsu(const struct quirc *q)
 			break;
 
 		sumB += i * histogram[i];
-		const double m1 = sumB / q1;
+		const quirc_float_t m1 = sumB / q1;
-		const double m2 = (sum - sumB) / q2;
+		const quirc_float_t m2 = (sum - sumB) / q2;
-		const double m1m2 = m1 - m2;
+		const quirc_float_t m1m2 = m1 - m2;
-		const double variance = m1m2 * m1m2 * q1 * q2;
+		const quirc_float_t variance = m1m2 * m1m2 * q1 * q2;
 		if (variance >= max) {
 			threshold = i;
 			max = variance;
@@ -582,7 +586,7 @@ static void find_alignment_pattern(struct quirc *q, int index)
 	int size_estimate;
 	int step_size = 1;
 	int dir = 0;
-	double u, v;
+	quirc_float_t u, v;
 
 	/* Grab our previous estimate of the alignment pattern corner */
 	memcpy(&b, &qr->align, sizeof(b));
@@ -648,10 +652,10 @@ static void find_leftmost_to_line(void *user_data, int y, int left, int right)
 	}
 }
 
-static double length(struct quirc_point a, struct quirc_point b)
+static quirc_float_t length(struct quirc_point a, struct quirc_point b)
 {
-	double x = abs(a.x - b.x) + 1;
+	quirc_float_t x = abs(a.x - b.x) + 1;
-	double y = abs(a.y - b.y) + 1;
+	quirc_float_t y = abs(a.y - b.y) + 1;
 	return sqrt(x * x +  y * y);
 }
 /* Estimate grid size by determing distance between capstones
@@ -664,15 +668,15 @@ static void measure_grid_size(struct quirc *q, int index)
 	struct quirc_capstone *b = &(q->capstones[qr->caps[1]]);
 	struct quirc_capstone *c = &(q->capstones[qr->caps[2]]);
 
-	double ab = length(b->corners[0], a->corners[3]);
+	quirc_float_t ab = length(b->corners[0], a->corners[3]);
-	double capstone_ab_size = (length(b->corners[0], b->corners[3]) + length(a->corners[0], a->corners[3]))/2.0;
+	quirc_float_t capstone_ab_size = (length(b->corners[0], b->corners[3]) + length(a->corners[0], a->corners[3]))/2.0;
-	double ver_grid = 7.0 * ab / capstone_ab_size;
+	quirc_float_t ver_grid = 7.0 * ab / capstone_ab_size;
 
-	double bc = length(b->corners[0], c->corners[1]);
+	quirc_float_t bc = length(b->corners[0], c->corners[1]);
-	double capstone_bc_size = (length(b->corners[0], b->corners[1]) + length(c->corners[0], c->corners[1]))/2.0;
+	quirc_float_t capstone_bc_size = (length(b->corners[0], b->corners[1]) + length(c->corners[0], c->corners[1]))/2.0;
-	double hor_grid = 7.0 * bc / capstone_bc_size;
+	quirc_float_t hor_grid = 7.0 * bc / capstone_bc_size;
 	
-	double grid_size_estimate = (ver_grid + hor_grid) / 2;
+	quirc_float_t grid_size_estimate = (ver_grid + hor_grid) / 2;
 
 	int ver = (int)((grid_size_estimate - 17.0 + 2.0) / 4.0);
 	
@@ -703,7 +707,7 @@ static int fitness_cell(const struct quirc *q, int index, int x, int y)
 
 	for (v = 0; v < 3; v++)
 		for (u = 0; u < 3; u++) {
-			static const double offsets[] = {0.3, 0.5, 0.7};
+			static const quirc_float_t offsets[] = {0.3, 0.5, 0.7};
 			struct quirc_point p;
 
 			perspective_map(qr->c, x + offsets[u],
@@ -806,7 +810,7 @@ static void jiggle_perspective(struct quirc *q, int index)
 	struct quirc_grid *qr = &q->grids[index];
 	int best = fitness_all(q, index);
 	int pass;
-	double adjustments[8];
+	quirc_float_t adjustments[8];
 	int i;
 
 	for (i = 0; i < 8; i++)
@@ -816,9 +820,9 @@ static void jiggle_perspective(struct quirc *q, int index)
 		for (i = 0; i < 16; i++) {
 			int j = i >> 1;
 			int test;
-			double old = qr->c[j];
+			quirc_float_t old = qr->c[j];
-			double step = adjustments[j];
+			quirc_float_t step = adjustments[j];
-			double new;
+			quirc_float_t new;
 
 			if (i & 1)
 				new = old + step;
@@ -998,7 +1002,7 @@ fail:
 
 struct neighbour {
 	int		index;
-	double		distance;
+	quirc_float_t		distance;
 };
 
 struct neighbour_list {
@@ -1015,7 +1019,7 @@ static void test_neighbours(struct quirc *q, int i,
 		const struct neighbour *hn = &hlist->n[j];
 		for (int k = 0; k < vlist->count; k++) {
 			const struct neighbour *vn = &vlist->n[k];
-			double squareness = fabs(1.0 - hn->distance / vn->distance);
+			quirc_float_t squareness = fabs(1.0 - hn->distance / vn->distance);
 			if (squareness < 0.2)
 				record_qr_grid(q, hn->index, i, vn->index);
 		}
@@ -1037,7 +1041,7 @@ static void test_grouping(struct quirc *q, unsigned int i)
 	 */
 	for (j = 0; j < q->num_capstones; j++) {
 		struct quirc_capstone *c2 = &q->capstones[j];
-		double u, v;
+		quirc_float_t u, v;
 
 		if (i == j)
 			continue;

lib/quirc_internal.h

@@ -46,6 +46,20 @@ typedef uint16_t quirc_pixel_t;
 #error "QUIRC_MAX_REGIONS > 65534 is not supported"
 #endif
 
+#ifdef QUIRC_FLOAT_TYPE
+/* Quirc uses double precision floating point internally by default.
+ * On platforms with a single precision FPU but no double precision FPU,
+ * this can be changed to float by defining QUIRC_FLOAT_TYPE.
+ *
+ * When setting QUIRC_FLOAT_TYPE to 'float', consider also defining QUIRC_USE_TGMATH.
+ * This will use the type-generic math functions (tgmath.h, C99 or later) instead of the normal ones,
+ * which will allow the compiler to use the correct overloaded functions for the type.
+ */
+typedef QUIRC_FLOAT_TYPE quirc_float_t;
+#else
+typedef double quirc_float_t;
+#endif
+
 struct quirc_region {
 	struct quirc_point	seed;
 	int			count;
@@ -58,7 +72,7 @@ struct quirc_capstone {
 
 	struct quirc_point	corners[4];
 	struct quirc_point	center;
-	double			c[QUIRC_PERSPECTIVE_PARAMS];
+	quirc_float_t		c[QUIRC_PERSPECTIVE_PARAMS];
 
 	int			qr_grid;
 };
@@ -76,7 +90,7 @@ struct quirc_grid {
 
 	/* Grid size and perspective transform */
 	int			grid_size;
-	double			c[QUIRC_PERSPECTIVE_PARAMS];
+	quirc_float_t		c[QUIRC_PERSPECTIVE_PARAMS];
 };
 
 struct quirc_flood_fill_vars {