commit 99be43e366b5c4d9dbaeaad9d885ff6d7f3aec4d
Author:     Mattias Andrée <maandree_AT_kth.se>
AuthorDate: Fri Jan 20 00:34:09 2017 +0100
Commit:     Mattias Andrée <maandree_AT_kth.se>
CommitDate: Fri Jan 20 01:13:46 2017 +0100
    Fix and improve blind-gauss-blur, and fix new bug in blind-from-image
    
    Signed-off-by: Mattias Andrée <maandree_AT_kth.se>
diff --git a/TODO b/TODO
index 3ff7bd7..a0d8f11 100644
--- a/TODO
+++ b/TODO
_AT_@ -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
index b66d127..83ed801 100644
--- a/src/blind-from-image.c
+++ b/src/blind-from-image.c
_AT_@ -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
index 8ba9743..a81a3ec 100644
--- a/src/blind-gauss-blur.c
+++ b/src/blind-gauss-blur.c
_AT_@ -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];
 
_AT_@ -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;
_AT_@ -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);
 
_AT_@ -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;
_AT_@ -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);
_AT_@ -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
index 2ad30ed..04dd516 100644
--- a/src/stream.c
+++ b/src/stream.c
_AT_@ -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
index 5fa3eea..f2d88ed 100644
--- a/src/stream.h
+++ b/src/stream.h
_AT_@ -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
index b87e9b5..c315010 100644
--- a/src/util.c
+++ b/src/util.c
_AT_@ -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
index bb9a8e4..d45433c 100644
--- a/src/util/jobs.h
+++ b/src/util/jobs.h
_AT_@ -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);
Received on Fri Jan 20 2017 - 01:14:26 CET
This archive was generated by hypermail 2.3.0
: Fri Jan 20 2017 - 01:24:20 CET