diff options
Diffstat (limited to 'src/blind-invert-matrix.c')
| -rw-r--r-- | src/blind-invert-matrix.c | 118 |
1 files changed, 61 insertions, 57 deletions
diff --git a/src/blind-invert-matrix.c b/src/blind-invert-matrix.c index 9befe67..fb86cc1 100644 --- a/src/blind-invert-matrix.c +++ b/src/blind-invert-matrix.c @@ -1,62 +1,20 @@ /* See LICENSE file for copyright and license details. */ +#ifndef TYPE #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++)\ - p2[c][0] -= p1[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;\ - 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 = 1 / 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); } +#define FILE "blind-invert-matrix.c" +#include "define-functions.h" int main(int argc, char *argv[]) { struct stream stream; - size_t width, x, y, row_size; - char *buf, *one = alloca(4 * sizeof(double)), *p; + size_t width, x, y, i, row_size; + char *buf, *one, *p; void (*process)(struct stream *stream, void *buf); ARGBEGIN { @@ -92,9 +50,8 @@ main(int argc, char *argv[]) eprintf("pixel format %s is not supported, try xyza\n", stream.pixfmt); } - memcpy(one + 1 * stream.chan_size, one, stream.chan_size); - memcpy(one + 2 * stream.chan_size, one, stream.chan_size); - memcpy(one + 3 * stream.chan_size, one, stream.chan_size); + for (i = 1; i < stream.n_chan; i++) + memcpy(one + i * stream.chan_size, one, stream.chan_size); width = stream.width > stream.height ? stream.width : 2 * stream.height; buf = emalloc2(width, stream.col_size); @@ -109,19 +66,17 @@ main(int argc, char *argv[]) } } if (equal) { - process(&stream, buf + 1 * stream.chan_size); + process(&stream, buf); 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, p + stream.chan_size, stream.chan_size); - memcpy(p + 2 * stream.chan_size, p, 2 * stream.chan_size); + for (i = 1; i < stream.n_chan; i++) + memcpy(p + i * stream.chan_size, p, stream.chan_size); } } } else { - process(&stream, buf + 0 * stream.chan_size); - process(&stream, buf + 1 * stream.chan_size); - process(&stream, buf + 2 * stream.chan_size); - process(&stream, buf + 3 * stream.chan_size); + for (i = 0; i < stream.n_chan; 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>"); @@ -130,3 +85,52 @@ main(int argc, char *argv[]) free(buf); return 0; } + +#else + +#define SUB_ROWS()\ + do {\ + p2 = matrix + r2 * cn;\ + t = p2[r1][0];\ + for (c = 0; c < cn; c++)\ + p2[c][0] -= p1[c][0] * t;\ + } while (0) + +static void +PROCESS(struct stream *stream, void *buf) +{ + 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; + 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 = 1 / 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(); + } +} + +#undef SUB_ROWS +#endif |