blind-colour-matrix.c (3441B)
1 /* See LICENSE file for copyright and license details. */ 2 #include "common.h" 3 4 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])") 5 6 static void 7 invert(double M[3][6]) 8 { 9 size_t r1, r2, i; 10 double t; 11 for (r1 = 0; r1 < 3; r1++) { 12 if (!M[r1][r1]) { 13 for (r2 = r1 + 1; r2 < 3 && !M[r2][r1]; r2++); 14 if (r2 >= 3) 15 eprintf("the colour space's rank is less than 3\n"); 16 for (i = 0; i < 6; i++) 17 t = M[r1][i], M[r1][i] = M[r2][i], M[r2][i] = t; 18 } 19 t = 1. / M[r1][r1]; 20 for (i = 0; i < 6; i++) 21 M[r1][i] *= t; 22 for (r2 = r1 + 1; r2 < 3; r2++) 23 for (i = 0, t = M[r2][r1]; i < 6; i++) 24 M[r2][i] -= M[r1][i] * t; 25 } 26 for (r1 = 3; --r1;) 27 for (r2 = r1; r2--;) 28 for (i = 0, t = M[r2][r1]; i < 6; i++) 29 M[r2][i] -= M[r1][i] * t; 30 } 31 32 int 33 main(int argc, char *argv[]) 34 { 35 static struct stream stream = { .width = 3, .height = 3, .frames = 1 }; 36 const char *pixfmt = "xyza"; 37 int ciexyy = 0; 38 double x[4], y[4], z[4], M[3][6], t; 39 double Mlf[9 * 4]; 40 float Mf[9 * 4]; 41 size_t i, j; 42 43 ARGBEGIN { 44 case 'F': 45 pixfmt = UARGF(); 46 break; 47 case 'z': 48 ciexyy = 1; 49 break; 50 default: 51 usage(); 52 } ARGEND; 53 54 if (argc != (3 - ciexyy) * 3 && argc != (3 - ciexyy) * 4) 55 usage(); 56 57 if (ciexyy) { 58 x[0] = etolf_arg("x1", argv[0]); 59 y[0] = etolf_arg("y1", argv[1]); 60 x[1] = etolf_arg("x2", argv[2]); 61 y[1] = etolf_arg("y2", argv[3]); 62 x[2] = etolf_arg("x3", argv[4]); 63 y[2] = etolf_arg("y3", argv[5]); 64 x[3] = argc > 6 ? etolf_arg("white-x", argv[6]) : D65_XYY_X; 65 y[3] = argc > 6 ? etolf_arg("white-y", argv[7]) : D65_XYY_Y; 66 for (i = 0; i < 4; i++) { 67 if (y[i]) { 68 z[i] = (1. - x[i] - y[i]) / y[i]; 69 x[i] /= y[i]; 70 y[i] = 1.; 71 } else { 72 x[i] = y[i] = z[i] = 1.; 73 } 74 } 75 } else { 76 x[0] = etolf_arg("X1", argv[0]); 77 y[0] = etolf_arg("Y1", argv[1]); 78 z[0] = etolf_arg("Z1", argv[2]); 79 x[1] = etolf_arg("X2", argv[3]); 80 y[1] = etolf_arg("Y2", argv[4]); 81 z[1] = etolf_arg("Z2", argv[5]); 82 x[2] = etolf_arg("X3", argv[6]); 83 y[2] = etolf_arg("Y3", argv[7]); 84 z[2] = etolf_arg("Z3", argv[8]); 85 x[3] = argc > 9 ? etolf_arg("white-X", argv[9]) : D65_XYZ_X; 86 y[3] = argc > 9 ? etolf_arg("white-Y", argv[10]) : 1; 87 z[3] = argc > 9 ? etolf_arg("white-Z", argv[11]) : D65_XYZ_Z; 88 } 89 90 for (i = 0; i < 3; i++) { 91 M[0][i] = x[i]; 92 M[1][i] = y[i]; 93 M[2][i] = z[i]; 94 M[i][3] = M[i][4] = M[i][5] = 0.; 95 M[i][3 + i] = 1.; 96 } 97 98 invert(M); 99 100 for (i = 0; i < 3; i++) { 101 t = M[i][3] * x[3] + M[i][4] * y[3] + M[i][5] * z[3]; 102 M[0][i] = t * x[i]; 103 M[1][i] = t * y[i]; 104 M[2][i] = t * z[i]; 105 } 106 107 for (i = 0; i < 3; i++) { 108 M[i][3] = M[i][4] = M[i][5] = 0.; 109 M[i][3 + i] = 1.; 110 } 111 112 invert(M); 113 114 eset_pixel_format(&stream, pixfmt); 115 fprint_stream_head(stdout, &stream); 116 efflush(stdout, "<stdout>"); 117 118 for (i = 0; i < 3; i++) { 119 for (j = 0; j < 3; j++) { 120 Mlf[i * 12 + j * 4 + 0] = M[i][3 + j]; 121 Mlf[i * 12 + j * 4 + 1] = M[i][3 + j]; 122 Mlf[i * 12 + j * 4 + 2] = M[i][3 + j]; 123 Mlf[i * 12 + j * 4 + 3] = 1.; 124 } 125 } 126 127 CHECK_ALPHA_CHAN(&stream); 128 CHECK_COLOUR_SPACE(&stream, CIEXYZ); 129 if (stream.encoding == DOUBLE) { 130 ewriteall(STDOUT_FILENO, Mlf, sizeof(Mlf), "<stdout>"); 131 } else if (stream.encoding == FLOAT) { 132 for (i = 0; i < ELEMENTSOF(Mlf); i++) 133 Mf[i] = (float)Mlf[i]; 134 ewriteall(STDOUT_FILENO, Mf, sizeof(Mf), "<stdout>"); 135 } else { 136 eprintf("pixel format %s is not supported, try xyza\n", stream.pixfmt); 137 } 138 139 return 0; 140 }