commit 99be43e366b5c4d9dbaeaad9d885ff6d7f3aec4d
parent e4880447680f6576fe571ac2c4bb0c3957b20ee9
Author: Mattias Andrée <maandree@kth.se>
Date: Fri, 20 Jan 2017 00:34:09 +0100
Fix and improve blind-gauss-blur, and fix new bug in blind-from-image
Signed-off-by: Mattias Andrée <maandree@kth.se>
Diffstat:
7 files changed, 309 insertions(+), 113 deletions(-)
diff --git a/TODO b/TODO
@@ -15,6 +15,5 @@ UNTESTED:
blind-cut
blind-extend
blind-from-text
- blind-gauss-blur
blind-rewrite-head
blind-to-text
diff --git a/src/blind-from-image.c b/src/blind-from-image.c
@@ -255,7 +255,7 @@ after_fork:
eprintf("%s\n", conv_fail_msg);
if (!headless) {
- FPRINTF_HEAD(stdout, 1, width, height, "xyza");
+ FPRINTF_HEAD_FMT(stdout, "%i", 1, "%s", width, "%s", height, "xyza");
efflush(stdout, "<stdout>");
}
diff --git a/src/blind-gauss-blur.c b/src/blind-gauss-blur.c
@@ -9,15 +9,20 @@
#include <string.h>
#include <unistd.h>
-USAGE("[-j jobs] [-ac] sd-stream")
+USAGE("[-j jobs] [-s spread | -s 'auto'] [-achvy] sd-stream")
static int chroma = 0;
static int noalpha = 0;
+static int vertical = 0;
+static int horizontal = 0;
+static int measure_y_only = 0;
+static int auto_spread = 0;
static size_t jobs = 1;
+static size_t spread = 0;
static void
process_xyza(char *restrict output, char *restrict cbuf, char *restrict sbuf,
- struct stream *colour, struct stream *sigma, size_t n, size_t sn)
+ struct stream *colour, struct stream *sigma, size_t cn, size_t sn)
{
typedef double pixel_t[4];
@@ -25,120 +30,264 @@ process_xyza(char *restrict output, char *restrict cbuf, char *restrict sbuf,
pixel_t *restrict sig = (pixel_t *)sbuf;
pixel_t *img = (pixel_t *)output;
pixel_t c, k;
- double *p;
size_t x1, y1, i1, x2, y2, i2;
- double dy, d, m, dblurred;
+ double d, m, X, Z;
int i, blurred, blur[3] = {0, 0, 0};
- size_t start = 0, end = colour->height;
+ size_t start, end, x2start, x2end, y2start, y2end;
int is_master;
+ pid_t *children;
- /* premultiply alpha channel */
- if (!noalpha) {
- for (y1 = i1 = 0; y1 < colour->height; y1++) {
- for (x1 = 0; x1 < colour->width; x1++, i1++) {
- clr[i1][0] *= clr[i1][3];
- clr[i1][1] *= clr[i1][3];
- clr[i1][2] *= clr[i1][3];
+ X = D65_XYY_X / D65_XYY_Y;
+ Z = 1 / D65_XYY_Y - 1 - X;
+
+ y2start = x2start = 0;
+ x2end = colour->width;
+ y2end = colour->height;
+
+ if (chroma || !noalpha) {
+ start = 0, end = colour->height;
+ is_master = efork_jobs(&start, &end, jobs, &children);
+
+ /* premultiply alpha channel */
+ if (!noalpha) {
+ i1 = start * colour->width;
+ for (y1 = start; y1 < end; y1++) {
+ for (x1 = 0; x1 < colour->width; x1++, i1++) {
+ clr[i1][0] *= clr[i1][3];
+ clr[i1][1] *= clr[i1][3];
+ clr[i1][2] *= clr[i1][3];
+ }
+ }
+ }
+
+ /* convert colour model */
+ if (chroma) {
+ i1 = start * colour->width;
+ for (y1 = start; y1 < end; y1++) {
+ for (x1 = 0; x1 < colour->width; x1++, i1++) {
+ clr[i1][0] = clr[i1][0] / X - clr[i1][1];
+ clr[i1][2] = clr[i1][2] / Z - clr[i1][1];
+ /*
+ * Explaination:
+ * Y is the luma and ((X / Xn - Y / Yn), (Z / Zn - Y / Yn))
+ * is the chroma (according to CIELAB), where (Xn, Yn, Zn)
+ * is the white point.
+ */
+ }
}
}
+ /* Conversion makes no difference if blur is applied to all
+ * parameters:
+ *
+ * Gaussian blur:
+ *
+ * ∞ ∞
+ * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²)
+ * V′ (x₀,y₀) = │ │ ────── e dxdy
+ * σ ⌡ ⌡ 2πσ²
+ * −∞ −∞
+ *
+ * With linear transformation, F:
+ *
+ * ∞ ∞
+ * ⌠ ⌠ F(V(x,y)) −((x−x₀)² + (y−y₀)²)/(2σ²)
+ * V′ (x₀,y₀) = F⁻¹ │ │ ───────── e dxdy
+ * σ ⌡ ⌡ 2πσ²
+ * −∞ −∞
+ *
+ * ∞ ∞
+ * ⌠ ⌠ ⎛V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²)⎞
+ * V′ (x₀,y₀) = F⁻¹ │ │ F⎜────── e ⎟ dxdy
+ * σ ⌡ ⌡ ⎝ 2πσ² ⎠
+ * −∞ −∞
+ *
+ * ∞ ∞
+ * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²)
+ * V′ (x₀,y₀) = (F⁻¹ ∘ F) │ │ ────── e dxdy
+ * σ ⌡ ⌡ 2πσ²
+ * −∞ −∞
+ *
+ * ∞ ∞
+ * ⌠ ⌠ V(x,y) −((x−x₀)² + (y−y₀)²)/(2σ²)
+ * V′ (x₀,y₀) = │ │ ────── e dxdy
+ * σ ⌡ ⌡ 2πσ²
+ * −∞ −∞
+ *
+ * Just like expected, the colour space should not affect the
+ * result of guassian blur as long as it is linear.
+ */
+
+ ejoin_jobs(is_master, children);
}
- is_master = efork_jobs(&start, &end, jobs);
+ /*
+ * This is not a regular simple gaussian blur implementation.
+ * This implementation is able to apply different levels of
+ * blur on different pixels. It's therefore written a bit
+ * oldly. Instead of going through each pixel and calculate
+ * the new value for each pixel, it goes through each pixel
+ * and smears it out to the other pixels.
+ */
+
+#define BLUR_PIXEL_PROLOGUE(DIR)\
+ if (sig[i1][3] == 0)\
+ goto no_blur_##DIR;\
+ if (chroma || measure_y_only) {\
+ k[0] = sig[i1][1] * sig[i1][3];\
+ if (auto_spread)\
+ spread = (size_t)(k[0] * 3 + 0.5);\
+ blur[2] = blur[0] = !!k[0];\
+ c[0] = k[0] *= k[0] * 2, c[0] = sqrt(c[0] * M_PI);\
+ k[0] = 1 / -k[0], c[0] = 1 / c[0];\
+ if (chroma) {\
+ k[2] = k[0];\
+ c[2] = c[0];\
+ c[1] = k[1] = 0;\
+ blur[1] = 0;\
+ } else {\
+ k[2] = k[1] = k[0];\
+ c[2] = c[1] = c[0];\
+ blur[1] = blur[0];\
+ }\
+ } else {\
+ if (auto_spread)\
+ spread = 0;\
+ for (i = 0; i < 3; i++) {\
+ k[i] = sig[i1][i] * sig[i1][3];\
+ if (auto_spread && spread < (size_t)(k[i] * 3 + 0.5))\
+ spread = (size_t)(k[i] * 3 + 0.5);\
+ blur[i] = !!k[i];\
+ c[i] = k[i] *= k[i] * 2, c[i] = sqrt(c[i] * M_PI);\
+ k[i] = 1 / -k[i], c[i] = 1 / c[i];\
+ }\
+ }\
+ if (blur[0] + blur[1] + blur[2] == 0)\
+ goto no_blur_##DIR;
+
+#define BLUR_PIXEL(START, LOOP, DISTANCE)\
+ if (k[0] == k[1] && k[1] == k[2]) {\
+ START;\
+ for (LOOP) {\
+ d = (DISTANCE);\
+ d *= d;\
+ m = c[0] * exp(d * k[0]);\
+ img[i2][0] += clr[i1][0] * m;\
+ img[i2][1] += clr[i1][1] * m;\
+ img[i2][2] += clr[i1][2] * m;\
+ if (!noalpha)\
+ img[i2][3] += clr[i1][3] * m;\
+ }\
+ } else {\
+ blurred = 0;\
+ for (i = 0; i < 3; i++) {\
+ if (blur[i])\
+ blurred += 1;\
+ else\
+ img[i1][i] += clr[i1][i];\
+ }\
+ for (i = 0; i < 3; i++) {\
+ if (!blur[i])\
+ continue;\
+ START;\
+ if (!noalpha) {\
+ for (LOOP) {\
+ d = (DISTANCE);\
+ d *= d;\
+ m = c[i] * exp(d * k[i]);\
+ img[i2][i] += clr[i1][i] * m;\
+ img[i2][3] += clr[i1][3] * m / blurred;\
+ }\
+ } else {\
+ for (LOOP) {\
+ d = (DISTANCE);\
+ d *= d;\
+ m = c[i] * exp(d * k[i]);\
+ img[i2][i] += clr[i1][i] * m;\
+ }\
+ }\
+ }\
+ }
+
+#define BLUR_PIXEL_EPILOGUE(DIR)\
+ continue;\
+ no_blur_##DIR:\
+ img[i1][0] = clr[i1][0];\
+ img[i1][1] = clr[i1][1];\
+ img[i1][2] = clr[i1][2];\
+ img[i1][3] = clr[i1][3];
+
+#define BLUR(DIR, SETSTART, SETEND, START, LOOP, DISTANCE)\
+ do {\
+ memset(img, 0, cn);\
+ start = 0, end = colour->height;\
+ is_master = efork_jobs(&start, &end, jobs, &children);\
+ i1 = start * colour->width;\
+ if (noalpha)\
+ for (y1 = start; y1 < end; y1++)\
+ for (x1 = 0; x1 < colour->width; x1++, i1++)\
+ img[i1][3] = 1;\
+ for (y1 = start; y1 < end; y1++) {\
+ for (x1 = 0; x1 < colour->width; x1++, i1++) {\
+ BLUR_PIXEL_PROLOGUE(DIR);\
+ if (spread) {\
+ SETSTART;\
+ SETEND;\
+ }\
+ BLUR_PIXEL(START, LOOP, DISTANCE);\
+ BLUR_PIXEL_EPILOGUE(DIR);\
+ }\
+ }\
+ ejoin_jobs(is_master, children);\
+ } while (0)
/* blur */
- i1 = start * colour->width;
- for (y1 = i1 = start; y1 < end; y1++) {
- for (x1 = 0; x1 < colour->width; x1++, i1++) {
- if (sig[i1][3] == 0)
- goto no_blur;
- if (!chroma) {
- for (i = 0; i < 3; i++) {
- k[i] = sig[i1][i] * sig[i1][3], c[i] = k[i] *= k[i] * 2, c[i] *= M_PI;
- k[i] = 1 / k[i], c[i] = -1 / c[i];
- blur[i] = !sig[i1][1];
+ if (horizontal)
+ BLUR(horizontal,
+ x2start = spread > x1 ? 0 : x1 - spread,
+ x2end = spread + 1 > colour->width - x1 ? colour->width : x1 + spread + 1,
+ i2 = y1 * colour->width + x2start,
+ x2 = x2start; x2 < x2end; (x2++, i2++),
+ (ssize_t)x1 - (ssize_t)x2);
+ if (horizontal && vertical)
+ memcpy(clr, img, cn);
+ if (vertical)
+ BLUR(vertical,
+ y2start = spread > y1 ? 0 : y1 - spread,
+ y2end = spread + 1 > colour->height - y1 ? colour->height : y1 + spread + 1,
+ i2 = y2start * colour->width + x1,
+ y2 = y2start; y2 < y2end; (y2++, i2 += colour->width),
+ (ssize_t)y1 - (ssize_t)y2);
+
+ if (chroma || !noalpha) {
+ start = 0, end = colour->height;
+ is_master = efork_jobs(&start, &end, jobs, &children);
+
+ /* convert back to CIE XYZ */
+ if (chroma) {
+ i1 = start * colour->width;
+ for (y1 = start; y1 < end; y1++) {
+ for (x1 = 0; x1 < colour->width; x1++, i1++) {
+ img[i1][0] = (img[i1][0] + img[i1][1]) * X;
+ img[i1][2] = (img[i1][2] + img[i1][1]) * Z;
}
- } else {
- k[1] = sig[i1][1] * sig[i1][3], c[1] = k[1] *= k[1] * 2, c[1] *= M_PI;
- k[1] = 1 / -k[1], c[1] = 1 / c[1];
- blur[1] = !sig[i1][1];
}
- if (blur[0] + blur[1] + blur[2] == 0)
- goto no_blur;
-
- p = img[i1];
- p[0] = p[1] = p[2] = 0;
- p[3] = noalpha;
- if (k[0] == k[1] && k[1] == k[2]) {
- for (y2 = i2 = 0; y2 < colour->height; y2++) {
- dy = (ssize_t)y1 - (ssize_t)y2;
- dy *= dy;
- for (x2 = 0; x2 < colour->width; x2++, i2++) {
- d = (ssize_t)x1 - (ssize_t)x2;
- d = d * d + dy;
- m = c[i1] * exp(d * k[i1]);
- for (i = noalpha ? 3 : 4; i--;)
- p[i] += clr[i2][i] * m;
- }
- }
- } else {
- blurred = 0;
- for (i = 0; i < n; i++) {
- if (!blur[i]) {
- p[i] = clr[i1][i];
+ }
+
+ /* unpremultiply alpha channel */
+ if (!noalpha) {
+ i1 = start * colour->width;
+ for (y1 = start; y1 < end; y1++) {
+ for (x1 = 0; x1 < colour->width; x1++, i1++) {
+ if (!img[i1][3])
continue;
- }
- for (y2 = i2 = 0; y2 < colour->height; y2++) {
- dy = (ssize_t)y1 - (ssize_t)y2;
- dy *= dy;
- if (!noalpha) {
- for (x2 = 0; x2 < colour->width; x2++, i2++) {
- d = (ssize_t)x1 - (ssize_t)x2;
- d = d * d + dy;
- m = c[i1] * exp(d * k[i1]);
- p[i] += clr[i2][i] * m;
- p[3] += clr[i2][3] * m;
- }
- } else {
- for (x2 = 0; x2 < colour->width; x2++, i2++) {
- d = (ssize_t)x1 - (ssize_t)x2;
- d = d * d + dy;
- m = c[i1] * exp(d * k[i1]);
- p[i] += clr[i2][i] * m;
- }
- }
- }
- blurred += 1;
- }
- if (!noalpha) {
- dblurred = blurred;
- for (y2 = i2 = 0; y2 < colour->height; y2++)
- for (x2 = 0; x2 < colour->width; x2++, i2++)
- p[3] /= dblurred;
+ img[i1][0] /= img[i1][3];
+ img[i1][1] /= img[i1][3];
+ img[i1][2] /= img[i1][3];
}
}
-
- continue;
- no_blur:
- img[i1][0] = clr[i1][0];
- img[i1][1] = clr[i1][1];
- img[i1][2] = clr[i1][2];
- img[i1][3] = clr[i1][3];
}
- }
- ejoin_jobs(is_master);
-
- /* unpremultiply alpha channel */
- if (!noalpha) {
- for (y1 = i1 = 0; y1 < colour->height; y1++) {
- for (x1 = 0; x1 < colour->width; x1++, i1++) {
- if (!img[i1][3])
- continue;
- img[i1][0] /= img[i1][3];
- img[i1][1] /= img[i1][3];
- img[i1][2] /= img[i1][3];
- }
- }
+ ejoin_jobs(is_master, children);
}
(void) sigma;
@@ -149,6 +298,7 @@ int
main(int argc, char *argv[])
{
struct stream colour, sigma;
+ char *arg;
void (*process)(char *restrict output, char *restrict cbuf, char *restrict sbuf,
struct stream *colour, struct stream *sigma, size_t cn, size_t sn);
@@ -159,9 +309,25 @@ main(int argc, char *argv[])
case 'c':
chroma = 1;
break;
+ case 'h':
+ horizontal = 1;
+ break;
+ case 'v':
+ vertical = 1;
+ break;
+ case 'y':
+ measure_y_only = 1;
+ break;
case 'j':
jobs = etozu_flag('j', EARG(), 1, SHRT_MAX);
break;
+ case 's':
+ arg = EARG();
+ if (!strcmp(arg, "auto"))
+ auto_spread = 1;
+ else
+ spread = etozu_flag('s', arg, 1, SIZE_MAX);
+ break;
default:
usage();
} ARGEND;
@@ -169,6 +335,9 @@ main(int argc, char *argv[])
if (argc != 1)
usage();
+ if (!vertical && !horizontal)
+ vertical = horizontal = 1;
+
colour.file = "<stdin>";
colour.fd = STDIN_FILENO;
einit_stream(&colour);
@@ -184,6 +353,11 @@ main(int argc, char *argv[])
echeck_compat(&colour, &sigma);
+ if (jobs > colour.height)
+ jobs = colour.height;
+
+ fprint_stream_head(stdout, &colour);
+ efflush(stdout, "<stdout>");
process_each_frame_two_streams(&colour, &sigma, STDOUT_FILENO, "<stdout>", process);
return 0;
diff --git a/src/stream.c b/src/stream.c
@@ -356,7 +356,5 @@ nprocess_each_frame_two_streams(int status, struct stream *left, struct stream *
enwriteall(status, output_fd, image, lframe_size, output_fname);
}
- free(lbuf);
- free(rbuf);
munmap(image, 2 * lframe_size + rframe_size);
}
diff --git a/src/stream.h b/src/stream.h
@@ -16,6 +16,10 @@
fprintf(FP, "%zu %zu %zu %s\n%cuivf",\
(size_t)(FRAMES), (size_t)(WIDTH), (size_t)(HEIGHT), PIXFMT, 0)
+#define FPRINTF_HEAD_FMT(FP, FFRAMES, FRAMES, FWIDTH, WIDTH, FHEIGHT, HEIGHT, PIXFMT)\
+ fprintf(FP, FFRAMES" "FWIDTH" "FHEIGHT" %s\n%cuivf",\
+ FRAMES, WIDTH, HEIGHT, PIXFMT, 0)
+
#define einit_stream(...) eninit_stream(1, __VA_ARGS__)
#define eset_pixel_size(...) enset_pixel_size(1, __VA_ARGS__)
#define eread_stream(...) enread_stream(1, __VA_ARGS__)
diff --git a/src/util.c b/src/util.c
@@ -167,35 +167,56 @@ enfork(int status)
}
+/* If <() is used in Bash (possibily other shells), that process becomes
+ * child of the process for each <() is used. Therefore, we cannot simply
+ * wait until the last child has been reaped, or even the expected number
+ * of children has been reaped, we must instead remember the PID of each
+ * child we created and wait for all of them to be reaped. { */
+
int
-enfork_jobs(int status, size_t *start, size_t *end, size_t jobs)
+enfork_jobs(int status, size_t *start, size_t *end, size_t jobs, pid_t **pids)
{
size_t j, s = *start, n = *end - *start;
- if (jobs < 2)
+ pid_t pid;
+ if (jobs < 2) {
+ *pids = NULL;
return 1;
+ }
*end = n / jobs + s;
+ *pids = enmalloc(status, jobs * sizeof(**pids));
for (j = 1; j < jobs; j++) {
- if (!enfork(status)) {
+ pid = enfork(status);
+ if (!pid) {
#if defined(HAVE_PRCTL) && defined(PR_SET_PDEATHSIG)
prctl(PR_SET_PDEATHSIG, SIGKILL);
#endif
*start = n * (j + 0) / jobs + s;
*end = n * (j + 1) / jobs + s;
return 0;
+ } else {
+ (*pids)[j - 1] = pid;
}
}
+ (*pids)[jobs - 1] = -1;
return 1;
}
void
-enjoin_jobs(int status, int is_master)
+enjoin_jobs(int status, int is_master, pid_t *pids)
{
int stat;
+ size_t i;
if (!is_master)
exit(0);
- while (wait(&stat) != -1)
- if (!stat)
+ if (!pids)
+ return;
+ for (i = 0; pids[i] != -1; i++) {
+ if (waitpid(pids[i], &stat, 0) == -1)
+ enprintf(status, "waitpid:");
+ if (stat)
exit(status);
- if (errno != ECHILD)
- enprintf(status, "wait:");
+ }
+ free(pids);
}
+
+/* } */
diff --git a/src/util/jobs.h b/src/util/jobs.h
@@ -3,5 +3,5 @@
#define efork_jobs(...) enfork_jobs(1, __VA_ARGS__)
#define ejoin_jobs(...) enjoin_jobs(1, __VA_ARGS__)
-int enfork_jobs(int status, size_t *start, size_t *end, size_t jobs);
-void enjoin_jobs(int status, int is_master);
+int enfork_jobs(int status, size_t *start, size_t *end, size_t jobs, pid_t **pids);
+void enjoin_jobs(int status, int is_master, pid_t *pids);
