Add gaussian blur and unshaping kernels

Signed-off-by: Mattias Andrée <maandree@kth.se>

2 files changed, 90 insertions, 5 deletions

diff --git a/src/blind-kernel.c b/src/blind-kernel.c
index 7ec22f8..f7405c1 100644
--- a/src/blind-kernel.c
+++ b/src/blind-kernel.c
@@ -2,18 +2,20 @@
 #include "stream.h"
 #include "util.h"
 
+#include <math.h>
 #include <string.h>
 
 USAGE("[-xyza] kernel [parameter] ...")
 
-#define SUBUSAGE(FORMAT) "Usage: %s [-xyza] " FORMAT, argv0
-#define MATRIX(...) ((const double[]){__VA_ARGS__});
+#define SUBUSAGE(FORMAT)      "Usage: %s [-xyza] " FORMAT, argv0
+#define MATRIX(...)           ((const double[]){__VA_ARGS__});
 #define STREQ3(A, B1, B2, B3) (!strcmp(A, B1) || !strcmp(A, B2) || !strcmp(A, B3))
 
 #define LIST_KERNELS\
-	X(kernel_kirsch, "kirsch")\
+	X(kernel_kirsch,   "kirsch")\
 	X(kernel_box_blur, "box blur")\
-	X(kernel_sharpen, "sharpen")
+	X(kernel_sharpen,  "sharpen")\
+	X(kernel_gaussian, "gaussian")
 
 static const double *
 kernel_kirsch(int argc, char *argv[], size_t *rows, size_t *cols, double **free_this)
@@ -55,7 +57,7 @@ kernel_box_blur(int argc, char *argv[], size_t *rows, size_t *cols, double **fre
 	}
 	*rows = 2 * sy + 1;
 	*cols = 2 * sx + 1;
-	*free_this = cells = emalloc2(*rows, *cols);
+	*free_this = cells = emalloc3(*rows, *cols, sizeof(double));
 	n = (2 * sy + 1) * (2 * sx + 1);
 	value = 1 / (double)n;
 	for (i = 0; i < n; i++)
@@ -74,6 +76,75 @@ kernel_sharpen(int argc, char *argv[], size_t *rows, size_t *cols, double **free
 	(void) argv;
 }
 
+static const double *
+kernel_gaussian(int argc, char *argv[], size_t *rows, size_t *cols, double **free_this)
+{
+	size_t spread = 0, y, x;
+	int unsharpen = 0;
+	double sigma, value, c, k;
+	char *arg;
+
+#define argv0 arg
+	argc++, argv--;
+	ARGBEGIN {
+	case 's':
+		if (!(arg = ARGF()))
+			goto usage;
+		spread = etozu_flag('s', arg, 1, SIZE_MAX / 2);
+		break;
+	case 'u':
+		unsharpen = 1;
+		break;
+	} ARGEND;
+#undef argv0
+
+	if (argc != 1)
+		goto usage;
+
+	sigma = etof_arg("standard-deviation", argv[0]);
+
+	if (!spread)
+		spread = (size_t)(sigma * 3.0 + 0.5);
+	*rows = *cols = spread * 2 + 1;
+
+	*free_this = emalloc3(*rows, *cols, sizeof(double));
+
+	k = sigma * sigma * 2;
+	c = M_PI * k;
+	c = sqrt(c);
+	c = 1.0 / c;
+	k = 1.0 / -k;
+
+	for (x = 0; x <= spread; x++) {
+		value = spread - x;
+		value *= value * k;
+		value = exp(value) * c;
+		for (y = 0; y < *rows; y++) {
+			(*free_this)[y * *cols + x] = value;
+			(*free_this)[y + 1 * *cols + *cols - 1 - x] = value;
+		}
+	}
+
+	for (y = 0; y <= spread; y++) {
+		value = spread - y;
+		value *= value * k;
+		value = exp(value) * c;
+		for (x = 0; x < *cols; x++) {
+			(*free_this)[y * *cols + x] *= value;
+			(*free_this)[y + 1 * *cols + *cols - 1 - x] *= value;
+		}
+	}
+
+	if (unsharpen)
+		(*free_this)[spread * *cols + spread] -= 2.0;
+
+	return *free_this;
+
+usage:
+	eprintf(SUBUSAGE("'gaussian' [-s spread] [-u] standard-deviation"));
+	return NULL;
+}
+
 int
 main(int argc, char *argv[])
 {
diff --git a/src/util/emalloc.h b/src/util/emalloc.h
index 214a773..39efcea 100644
--- a/src/util/emalloc.h
+++ b/src/util/emalloc.h
@@ -4,12 +4,14 @@
 
 #define emalloc(...)   enmalloc(1, __VA_ARGS__)
 #define emalloc2(...)  enmalloc2(1, __VA_ARGS__)
+#define emalloc3(...)  enmalloc3(1, __VA_ARGS__)
 #define ecalloc(...)   encalloc(1, __VA_ARGS__)
 #define erealloc(...)  enrealloc(1, __VA_ARGS__)
 #define erealloc2(...) enrealloc2(1, __VA_ARGS__)
 #define erealloc3(...) enrealloc3(1, __VA_ARGS__)
 
 #define malloc2(n, m)           malloc(n * m);
+#define malloc3(n1, n2, n3)     malloc(n1 * n2 * n3);
 #define realloc2(p, n, m)       realloc(p, n * m);
 #define realloc3(p, n1, n2, n3) realloc(p, n1 * n2 * n3);
 
@@ -32,6 +34,18 @@ enmalloc2(int status, size_t n, size_t m)
 }
 
 static inline void *
+enmalloc3(int status, size_t n1, size_t n2, size_t n3)
+{
+	void *ptr;
+	size_t n = n1;
+	if (n2 > SIZE_MAX / n ||
+	    n3 > SIZE_MAX / (n *= n2) ||
+	    !(ptr = malloc(n * n3)))
+		enprintf(status, "malloc: out of memory\n");
+	return ptr;
+}
+
+static inline void *
 encalloc(int status, size_t n, size_t m)
 {
 	void *ptr = calloc(n, m);