commit 4621e6f35cda10654d2bbfb4df579ec89fb7f99d
parent fc3bad11b3628c3b77df97327f1cdaf19d96aa9e
Author: Mattias Andrée <maandree@kth.se>
Date:   Fri, 14 Jul 2017 19:25:21 +0200
Add blind-invert-matrix
Signed-off-by: Mattias Andrée <maandree@kth.se>
Diffstat:
6 files changed, 175 insertions(+), 0 deletions(-)
diff --git a/Makefile b/Makefile
@@ -39,6 +39,7 @@ BIN =\
 	blind-hexagon-tessellation\
 	blind-interleave\
 	blind-invert-luma\
+	blind-invert-matrix\
 	blind-linear-gradient\
 	blind-make-kernel\
 	blind-matrix-orthoproject\
diff --git a/README b/README
@@ -111,6 +111,9 @@ UTILITIES
        blind-invert-luma(1)
               Invert the luminosity of a video
 
+       blind-invert-matrix(1)
+              Invert matrix-video
+
        blind-linear-gradient(1)
               Generate a video with a linear gradient
 
diff --git a/man/blind-invert-matrix.1 b/man/blind-invert-matrix.1
@@ -0,0 +1,35 @@
+.TH BLIND-INVERT-MATRIX 1 blind
+.SH NAME
+blind-invert-matrix - Invert matrix-vidoes
+.SH SYNOPSIS
+.B blind-invert-matrix
+[-e]
+.SH DESCRIPTION
+.B blind-invert-matrix
+reads a video representing a matrix from
+stdin, inverts it, and prints the results
+to stdout.
+.P
+The matrix must be at least as wide as it
+is tall. If the matrix is wider than it is
+tall, it is treated as an augmented matrix,
+and the unaugmented square matrix is
+eliminated to the identity matrix and the
+resuling augment is printed.
+.SH OPTIONS
+.TP
+.B -e
+Apply optimisation that assumes all channels
+are identical.
+.SH SEE ALSO
+.BR blind (7),
+.BR blind-multiply-matrices (1),
+.BR blind-transpose (1),
+.BR blind-flip (1),
+.BR blind-flop (1),
+.BR blind-rotate-90 (1),
+.BR blind-rotate-180 (1),
+.BR blind-rotate-270 (1)
+.SH AUTHORS
+Mattias Andrée
+.RI < maandree@kth.se >
diff --git a/man/blind-multiply-matrices.1 b/man/blind-multiply-matrices.1
@@ -38,6 +38,7 @@ in reverse order.
 .BR blind-matrix-shear (1),
 .BR blind-matrix-translate (1),
 .BR blind-matrix-transpose (1),
+.BR blind-invert-matrix (1),
 .BR blind-transpose (1),
 .BR blind-flip (1),
 .BR blind-flop (1),
diff --git a/man/blind.7 b/man/blind.7
@@ -127,6 +127,9 @@ Framewise interleave videos
 .BR blind-invert-luma (1)
 Invert the luminosity of a video
 .TP
+.BR blind-invert-matrix (1)
+Invert matrix-video
+.TP
 .BR blind-linear-gradient (1)
 Generate a video with a linear gradient
 .TP
diff --git a/src/blind-invert-matrix.c b/src/blind-invert-matrix.c
@@ -0,0 +1,132 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+USAGE("")
+
+static int equal = 0;
+
+#define SUB_ROWS()\
+	do {\
+		p2 = matrix + r2 * cn;\
+		t = p2[r1][0];\
+		for (c = 0; c < cn; c++)\
+			p1[c][0] -= p2[c][0] * t;\
+	} while (0)
+
+#define PROCESS(TYPE)\
+	do {\
+		typedef TYPE pixel_t[4];\
+		size_t rn = stream->height, r1, r2, c;\
+		size_t cn = stream->width > rn ? stream->width : 2 * rn;\
+		pixel_t *matrix = buf, *p1, *p2 = NULL;\
+		TYPE t;\
+		\
+		for (r1 = 0; r1 < rn; r1++) {\
+			p1 = matrix + r1 * cn;\
+ 			if (!p1[r1][0]) {\
+				for (r2 = r1 + 1; r2 < rn; r2++) {\
+					p2 = matrix + r2 * cn;\
+					if (p2[r1][0])\
+						break;\
+				}\
+				if (r2 == rn)\
+					eprintf("matrix is not invertable\n");\
+				for (c = 0; c < cn; c++)\
+					t = p1[c][0], p1[c][0] = p2[c][0], p2[c][0] = t;\
+			}\
+			t = p1[r1][0];\
+			for (c = 0; c < cn; c++)\
+				p1[c][0] /= t;\
+			for (r2 = r1 + 1; r2 < rn; r2++)\
+				SUB_ROWS();\
+		}\
+		\
+		for (r1 = rn; r1--;) {\
+			p1 = matrix + r1 * cn;\
+			for (r2 = r1; r2--;)\
+				SUB_ROWS();\
+		}\
+	} while (0)
+
+static void process_lf(struct stream *stream, void *buf) { PROCESS(double); }
+static void process_f (struct stream *stream, void *buf) { PROCESS(float); }
+
+int
+main(int argc, char *argv[])
+{
+	struct stream stream;
+	size_t width, x, y, row_size, chan_size;
+	char *buf, *one = alloca(4 * sizeof(double)), *p;
+	void (*process)(struct stream *stream, void *buf);
+
+	ARGBEGIN {
+	case 'e':
+		equal = 1;
+		break;
+	default:
+		usage();
+	} ARGEND;
+
+	if (argc)
+		usage();
+
+	eopen_stream(&stream, NULL);
+	echeck_dimensions(&stream, WIDTH | HEIGHT, NULL);
+	width = stream.width;
+	if (stream.width < stream.height)
+		eprintf("<stdin>: the video must be at least as wide as it is tall\n");
+	else if (stream.width > stream.height)
+		stream.width -= stream.height;
+	fprint_stream_head(stdout, &stream);
+	stream.width = width;
+	efflush(stdout, "<stdout>");
+
+	if (!strcmp(stream.pixfmt, "xyza")) {
+		*(double *)one = 1;
+		process = process_lf;
+	} else if (!strcmp(stream.pixfmt, "xyza f")) {
+		*(float *)one = 1;
+		process = process_f;
+	} else {
+		eprintf("pixel format %s is not supported, try xyza\n", stream.pixfmt);
+	}
+
+	chan_size = stream.pixel_size / 4;
+	memcpy(one + 1 * chan_size, one, chan_size);
+	memcpy(one + 2 * chan_size, one, chan_size);
+	memcpy(one + 3 * chan_size, one, chan_size);
+
+	width = stream.width > stream.height ? stream.width : 2 * stream.height;
+	buf = emalloc2(width, stream.col_size);
+	row_size = width * stream.pixel_size;
+
+	while (eread_frame(&stream, buf)) {
+		if (stream.width == stream.height) {
+			for (y = stream.height; y--;) {
+				memmove(buf + y * row_size, buf + y * stream.row_size, stream.row_size);
+				memset(buf + y * row_size + stream.row_size, 0, stream.row_size);
+				memcpy(buf + y * row_size + y * stream.pixel_size, one, stream.pixel_size);
+			}
+		}
+		if (equal) {
+			process(&stream, buf + 1 * 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;
+					memcpy(p + chan_size, p, chan_size);
+					memcpy(p + 2 * chan_size, p, 2 * chan_size);
+				}
+			}
+		} else {
+			process(&stream, buf + 0 * chan_size);
+			process(&stream, buf + 1 * chan_size);
+			process(&stream, buf + 2 * chan_size);
+			process(&stream, buf + 3 * chan_size);
+		}
+		for (y = 0; y < stream.height; y++)
+			ewriteall(STDOUT_FILENO, buf + y * row_size + stream.row_size, stream.row_size, "<stdout>");
+	}
+
+	free(buf);
+	return 0;
+}