Add blind-colour-matrix, and blind-invert-matrix: fix -e and add -axyz

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

2 files changed, 160 insertions, 8 deletions

diff --git a/src/blind-colour-matrix.c b/src/blind-colour-matrix.c
new file mode 100644
index 0000000..e4ff869
--- /dev/null
+++ b/src/blind-colour-matrix.c
@@ -0,0 +1,130 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+USAGE("[-F pixel-format] (-z x1 y1 x2 y2 x3 y3 white-x white-y | X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 [white-X white-Y white-Z])")
+
+int
+main(int argc, char *argv[])
+{
+	static struct stream stream = { .width = 3, .height = 3, .frames = 1 };
+	const char *pixfmt = "xyza";
+	int ciexyy = 0;
+	double x[4], y[4], z[4], M[3][6], t;
+	double Mlf[9 * 4];
+	float Mf[9 * 4];
+	size_t i, j, r1, r2;
+
+	ARGBEGIN {
+	case 'F':
+		pixfmt = UARGF();
+		break;
+	case 'z':
+		ciexyy = 1;
+		break;
+	default:
+		usage();
+	} ARGEND;
+
+	if (argc != (3 - ciexyy) * 3 && argc != (3 - ciexyy) * 4)
+		usage();
+
+	if (ciexyy) {
+		x[0] = etolf_arg("x1", argv[0]);
+		y[0] = etolf_arg("y1", argv[1]);
+		x[1] = etolf_arg("x2", argv[2]);
+		y[1] = etolf_arg("y2", argv[3]);
+		x[2] = etolf_arg("x3", argv[4]);
+		y[2] = etolf_arg("y3", argv[5]);
+		x[3] = argc > 6 ? etolf_arg("white-x", argv[6]) : D65_XYY_X;
+		y[3] = argc > 6 ? etolf_arg("white-y", argv[7]) : D65_XYY_Y;
+		for (i = 0; i < 4; i++) {
+			if (y[i]) {
+				z[i] = (1. - x[i] - y[i]) / y[i];
+				x[i] /= y[i];
+				y[i] = 1.;
+			} else {
+				x[i] = y[i] = z[i] = 1.;
+			}
+		}
+	} else {
+		x[0] = etolf_arg("X1", argv[0]);
+		y[0] = etolf_arg("Y1", argv[1]);
+		z[0] = etolf_arg("Z1", argv[2]);
+		x[1] = etolf_arg("X2", argv[3]);
+		y[1] = etolf_arg("Y2", argv[4]);
+		z[1] = etolf_arg("Z2", argv[5]);
+		x[2] = etolf_arg("X3", argv[6]);
+		y[2] = etolf_arg("Y3", argv[7]);
+		z[2] = etolf_arg("Z3", argv[8]);
+		x[3] = argc > 9 ? etolf_arg("white-X", argv[9])  : D65_XYZ_X;
+		y[3] = argc > 9 ? etolf_arg("white-Y", argv[10]) : 1;
+		z[3] = argc > 9 ? etolf_arg("white-Z", argv[11]) : D65_XYZ_Z;
+		for (i = 0; i < 4; i++) {
+			if (y[i] && y[i] != 1.) {
+				x[i] /= y[i];
+				z[i] /= y[i];
+				y[i] = 1.;
+			} else if (!y[i]) {
+				x[i] = y[i] = z[i] = 0.;
+			}
+		}
+	}
+
+	for (i = 0; i < 3; i++) {
+		M[0][i] = x[i];
+		M[1][i] = y[i];
+		M[2][i] = z[i];
+		M[i][3] = M[i][4] = M[i][5] = 0.;
+		M[i][3 + i] = 1.;
+	}
+
+	for (r1 = 0; r1 < 3; r1++) {
+		if (!M[r1][r1]) {
+			for (r2 = r1 + 1; r2 < 3 && !M[r2][r1]; r2++);
+			if (r2 >= 3)
+				eprintf("the colour space's rank is less than 3\n");
+			for (i = 0; i < 6; i++)
+				t = M[r1][i], M[r1][i] = M[r2][i], M[r2][i] = t;
+		}
+		t = 1. / M[r1][r1];
+		for (i = 0; i < 6; i++)
+			M[r1][i] *= t;
+		for (r2 = r1; r2--;)
+			for (i = 0, t = M[r2][r1]; i < 6; i++)
+				M[r2][i] -= M[r1][i] * t;
+	}
+	for (r1 = 3; r1--;)
+		for (r2 = r1; r2--;)
+			for (i = 0, t = M[r2][r1]; i < 6; i++)
+				M[r2][i] -= M[r1][i] * t;
+
+	for (i = 0; i < 3; i++) {
+		t = M[i][3] * x[3] + M[i][4] * y[3] + M[i][5] * z[3];
+		M[0][i] = t * x[i];
+		M[1][i] = t * y[i];
+		M[2][i] = t * z[i];
+	}
+
+	eset_pixel_format(&stream, pixfmt);
+	fprint_stream_head(stdout, &stream);
+	efflush(stdout, "<stdout>");
+
+	for (i = 0; i < 3; i++) {
+		for (j = 0; j < 3; j++) {
+			Mlf[i * 12 + j * 4 + 0] = M[i][j];
+			Mlf[i * 12 + j * 4 + 1] = M[i][j];
+			Mlf[i * 12 + j * 4 + 2] = M[i][j];
+			Mlf[i * 12 + j * 4 + 3] = 1.;
+		}
+	}
+
+	if (stream.encoding == DOUBLE) {
+		ewriteall(STDOUT_FILENO, Mlf, sizeof(Mlf), "<stdout>");
+	} else {
+		for (i = 0; i < ELEMENTSOF(Mlf); i++)
+			Mf[i] = (float)Mlf[i];
+		ewriteall(STDOUT_FILENO, Mf, sizeof(Mf), "<stdout>");
+	}
+
+	return 0;
+}
diff --git a/src/blind-invert-matrix.c b/src/blind-invert-matrix.c
index fb86cc1..66b1176 100644
--- a/src/blind-invert-matrix.c
+++ b/src/blind-invert-matrix.c
@@ -2,9 +2,11 @@
 #ifndef TYPE
 #include "common.h"
 
-USAGE("")
+USAGE("[-aexyz]")
 
 static int equal = 0;
+static int skip_ch[] = {0, 0, 0, 0};
+static size_t first_included = 0;
 
 #define FILE "blind-invert-matrix.c"
 #include "define-functions.h"
@@ -14,13 +16,25 @@ main(int argc, char *argv[])
 {
 	struct stream stream;
 	size_t width, x, y, i, row_size;
-	char *buf, *one, *p;
+	char *buf, *one, *p, *q;
 	void (*process)(struct stream *stream, void *buf);
 
 	ARGBEGIN {
+	case 'a':
+		skip_ch[3] = 1;
+		break;
 	case 'e':
 		equal = 1;
 		break;
+	case 'x':
+		skip_ch[0] = 1;
+		break;
+	case 'y':
+		skip_ch[1] = 1;
+		break;
+	case 'z':
+		skip_ch[2] = 1;
+		break;
 	default:
 		usage();
 	} ARGEND;
@@ -28,6 +42,11 @@ main(int argc, char *argv[])
 	if (argc)
 		usage();
 
+	while (first_included < ELEMENTSOF(skip_ch) && skip_ch[first_included])
+		first_included++;
+	if (first_included == ELEMENTSOF(skip_ch))
+		equal = 0;
+
 	eopen_stream(&stream, NULL);
 	echeck_dimensions(&stream, WIDTH | HEIGHT, NULL);
 	width = stream.width;
@@ -66,17 +85,20 @@ main(int argc, char *argv[])
 			}
 		}
 		if (equal) {
-			process(&stream, buf);
+			process(&stream, buf + first_included * stream.chan_size);
 			for (y = 0; y < stream.height; y++) {
 				for (x = 0; x < stream.width; x++) {
-					p = buf + y * row_size + x * stream.pixel_size;
-					for (i = 1; i < stream.n_chan; i++)
-						memcpy(p + i * stream.chan_size, p, stream.chan_size);
+					p = buf + y * row_size + x * stream.pixel_size + stream.col_size;
+					q = p + first_included * stream.chan_size;
+					for (i = 0; i < stream.n_chan; i++, p += stream.chan_size)
+						if (i != first_included && !skip_ch[i])
+							memcpy(p, q, stream.chan_size);
 				}
 			}
 		} else {
 			for (i = 0; i < stream.n_chan; i++)
-				process(&stream, buf + i * stream.chan_size);
+				if (!skip_ch[i])
+					process(&stream, buf + i * stream.chan_size);
 		}
 		for (y = 0; y < stream.height; y++)
 			ewriteall(STDOUT_FILENO, buf + y * row_size + stream.col_size, row_size - stream.col_size, "<stdout>");
@@ -107,7 +129,7 @@ PROCESS(struct stream *stream, void *buf)
 
 	for (r1 = 0; r1 < rn; r1++) {
 		p1 = matrix + r1 * cn;
-			if (!p1[r1][0]) {
+		if (!p1[r1][0]) {
 			for (r2 = r1 + 1; r2 < rn; r2++) {
 				p2 = matrix + r2 * cn;
 				if (p2[r1][0])