[hackers] [blind] Add gaussian blur and unshaping kernels || Mattias Andrée

[git]

commit 9e524b9dedeeea7f5d5e39edb80c588b161deed4
Author:     Mattias Andrée <maand...@kth.se>
AuthorDate: Sun May 14 00:23:52 2017 +0200
Commit:     Mattias Andrée <maand...@kth.se>
CommitDate: Sun May 14 00:23:52 2017 +0200

    Add gaussian blur and unshaping kernels
    
    Signed-off-by: Mattias Andrée <maand...@kth.se>

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);