ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/macbethcal.c
Revision: 2.7
Committed: Fri Oct 20 15:26:11 1995 UTC (28 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +223 -68 lines
Log Message:
added -c option for color input and changed option ordering

File Contents

# Content
1 /* Copyright (c) 1995 Regents of the University of California */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ LBL";
5 #endif
6
7 /*
8 * Calibrate a scanned MacBeth Color Checker Chart
9 *
10 * Produce a .cal file suitable for use with pcomb.
11 */
12
13 #include <stdio.h>
14 #ifdef MSDOS
15 #include <fcntl.h>
16 #endif
17 #include "color.h"
18 #include "resolu.h"
19 #include "pmap.h"
20
21 /* MacBeth colors */
22 #define DarkSkin 0
23 #define LightSkin 1
24 #define BlueSky 2
25 #define Foliage 3
26 #define BlueFlower 4
27 #define BluishGreen 5
28 #define Orange 6
29 #define PurplishBlue 7
30 #define ModerateRed 8
31 #define Purple 9
32 #define YellowGreen 10
33 #define OrangeYellow 11
34 #define Blue 12
35 #define Green 13
36 #define Red 14
37 #define Yellow 15
38 #define Magenta 16
39 #define Cyan 17
40 #define White 18
41 #define Neutral8 19
42 #define Neutral65 20
43 #define Neutral5 21
44 #define Neutral35 22
45 #define Black 23
46 /* computed from 5nm spectral measurements */
47 /* CIE 1931 2 degree obs, equal-energy white */
48 float mbxyY[24][3] = {
49 {0.462, 0.3769, 0.0932961}, /* DarkSkin */
50 {0.4108, 0.3542, 0.410348}, /* LightSkin */
51 {0.2626, 0.267, 0.181554}, /* BlueSky */
52 {0.36, 0.4689, 0.108447}, /* Foliage */
53 {0.2977, 0.2602, 0.248407}, /* BlueFlower */
54 {0.2719, 0.3485, 0.401156}, /* BluishGreen */
55 {0.52, 0.4197, 0.357899}, /* Orange */
56 {0.229, 0.1866, 0.103911}, /* PurplishBlue */
57 {0.4909, 0.3262, 0.242615}, /* ModerateRed */
58 {0.3361, 0.2249, 0.0600102}, /* Purple */
59 {0.3855, 0.4874, 0.42963}, /* YellowGreen */
60 {0.4853, 0.4457, 0.476343}, /* OrangeYellow */
61 {0.2026, 0.1369, 0.0529249}, /* Blue */
62 {0.3007, 0.4822, 0.221226}, /* Green */
63 {0.5805, 0.3238, 0.162167}, /* Red */
64 {0.4617, 0.472, 0.64909}, /* Yellow */
65 {0.4178, 0.2625, 0.233662}, /* Magenta */
66 {0.2038, 0.2508, 0.167275}, /* Cyan */
67 {0.3358, 0.337, 0.916877}, /* White */
68 {0.3338, 0.3348, 0.604678}, /* Neutral.8 */
69 {0.3333, 0.3349, 0.364566}, /* Neutral.65 */
70 {0.3353, 0.3359, 0.200238}, /* Neutral.5 */
71 {0.3363, 0.336, 0.0878721}, /* Neutral.35 */
72 {0.3346, 0.3349, 0.0308383} /* Black */
73 };
74
75 COLOR mbRGB[24]; /* MacBeth RGB values */
76
77 #define NMBNEU 6 /* Number of MacBeth neutral colors */
78 short mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};
79
80 #define NMBMOD 16 /* Number of MacBeth unsaturated colors */
81 short mbmod[NMBMOD] = {
82 DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen,
83 PurplishBlue,ModerateRed,YellowGreen,OrangeYellow,
84 Black,Neutral35,Neutral5,Neutral65,Neutral8,White
85 };
86
87 #define NMBSAT 8 /* Number of MacBeth saturated colors */
88 short mbsat[NMBSAT] = {
89 Red,Green,Blue,Magenta,Yellow,Cyan,
90 Orange,Purple
91 };
92
93 #define REQFLGS (1L<<White|1L<<Neutral8|1L<<Neutral65| \
94 1L<<Neutral5|1L<<Neutral35|1L<<Black)
95
96 #define CENTCVG 0.3 /* measured coverage of square sample */
97 #define FULLCVG 0.9 /* coverage of entire square */
98
99 int xmax, ymax; /* input image dimensions */
100 int bounds[4][2]; /* image coordinates of chart corners */
101 double imgxfm[3][3]; /* coordinate transformation matrix */
102
103 COLOR inpRGB[24]; /* measured or scanned input colors */
104 long inpflags = 0; /* flags of which colors were input */
105
106 COLOR bramp[NMBNEU][2]; /* brightness ramp (per primary) */
107 double solmat[3][3]; /* color mapping matrix */
108
109 FILE *debugfp = NULL; /* debug output picture */
110 char *progname;
111
112 extern char *malloc();
113
114
115 main(argc, argv)
116 int argc;
117 char **argv;
118 {
119 int inpispic = 1;
120 int i;
121
122 progname = argv[0];
123 for (i = 1; i < argc && argv[i][0] == '-'; i++)
124 switch (argv[i][1]) {
125 case 'd': /* debug output */
126 i++;
127 if (badarg(argc-i, argv+i, "s"))
128 goto userr;
129 if ((debugfp = fopen(argv[i], "w")) == NULL) {
130 perror(argv[i]);
131 exit(1);
132 }
133 #ifdef MSDOS
134 setmode(fileno(debugfp), O_BINARY);
135 #endif
136 newheader("RADIANCE", debugfp); /* start */
137 printargs(argc, argv, debugfp); /* header */
138 break;
139 case 'p': /* picture position */
140 if (badarg(argc-i-1, argv+i+1, "iiiiiiii"))
141 goto userr;
142 bounds[0][0] = atoi(argv[++i]);
143 bounds[0][1] = atoi(argv[++i]);
144 bounds[1][0] = atoi(argv[++i]);
145 bounds[1][1] = atoi(argv[++i]);
146 bounds[2][0] = atoi(argv[++i]);
147 bounds[2][1] = atoi(argv[++i]);
148 bounds[3][0] = atoi(argv[++i]);
149 bounds[3][1] = atoi(argv[++i]);
150 inpispic = 2;
151 break;
152 case 'c': /* color input */
153 inpispic = 0;
154 break;
155 default:
156 goto userr;
157 }
158 /* open files */
159 if (i < argc && freopen(argv[i], "r", stdin) == NULL) {
160 perror(argv[1]);
161 exit(1);
162 }
163 if (i+1 < argc && freopen(argv[i+1], "w", stdout) == NULL) {
164 perror(argv[2]);
165 exit(1);
166 }
167 if (inpispic) { /* load input picture header */
168 #ifdef MSDOS
169 setmode(fileno(stdin), O_BINARY);
170 #endif
171 if (checkheader(stdin, COLRFMT, NULL) < 0 ||
172 fgetresolu(&xmax, &ymax, stdin) < 0) {
173 fprintf(stderr, "%s: bad input picture\n", progname);
174 exit(1);
175 }
176 } else { /* else set default xmax and ymax */
177 xmax = 512;
178 ymax = 2*512/3;
179 }
180 if (inpispic != 2) { /* use default boundaries */
181 bounds[0][0] = bounds[2][0] = .029*xmax + .5;
182 bounds[0][1] = bounds[1][1] = .956*ymax + .5;
183 bounds[1][0] = bounds[3][0] = .971*xmax + .5;
184 bounds[2][1] = bounds[3][1] = .056*ymax + .5;
185 }
186 init(); /* initialize */
187 if (inpispic) /* get picture colors */
188 getpicture();
189 else
190 getcolors();
191 compute(); /* compute color mapping */
192 /* print comment */
193 printf("{ Color correction file computed by:\n\t");
194 printargs(argc, argv, stdout);
195 printf("}\n");
196 putmapping(); /* put out color mapping */
197 if (debugfp != NULL) /* put out debug picture */
198 if (inpispic)
199 picdebug();
200 else
201 clrdebug();
202 exit(0);
203 userr:
204 fprintf(stderr,
205 "Usage: %s [-d dbg.pic][-p xul yul xur yur xll yll xlr ylr] input.pic [output.cal]\n",
206 progname);
207 fprintf(stderr, " or: %s [-d dbg.pic] -c [xyY.dat [output.cal]]\n",
208 progname);
209 exit(1);
210 }
211
212
213 init() /* initialize */
214 {
215 double quad[4][2];
216 register int i;
217 /* make coordinate transformation */
218 quad[0][0] = bounds[0][0];
219 quad[0][1] = bounds[0][1];
220 quad[1][0] = bounds[1][0];
221 quad[1][1] = bounds[1][1];
222 quad[2][0] = bounds[3][0];
223 quad[2][1] = bounds[3][1];
224 quad[3][0] = bounds[2][0];
225 quad[3][1] = bounds[2][1];
226
227 if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
228 fprintf(stderr, "%s: bad chart boundaries\n", progname);
229 exit(1);
230 }
231 /* map MacBeth colors to RGB space */
232 for (i = 0; i < 24; i++)
233 xyY2RGB(mbRGB[i], mbxyY[i]);
234 }
235
236
237 int
238 chartndx(x, y, cvg) /* find color number for position */
239 int x, y;
240 double cvg;
241 {
242 double ipos[3], cpos[3];
243 int ix, iy;
244 double fx, fy;
245 double cmin, cmax;
246
247 ipos[0] = x;
248 ipos[1] = y;
249 ipos[2] = 1;
250 mx3d_transform(ipos, imgxfm, cpos);
251 cpos[0] /= cpos[2];
252 cpos[1] /= cpos[2];
253 if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
254 return(-1);
255 ix = cpos[0];
256 iy = cpos[1];
257 fx = cpos[0] - ix;
258 fy = cpos[1] - iy;
259 cmin = .5*(1.-cvg);
260 cmax = 1. - cmin;
261 if (fx < cmin || fx >= cmax || fy < cmin || fy >= cmax)
262 return(-1);
263 return(iy*6 + ix);
264 }
265
266
267 getpicture() /* load in picture colors */
268 {
269 COLR *scanln;
270 COLOR pval;
271 int ccount[24];
272 double d;
273 int y;
274 register int x, i;
275
276 scanln = (COLR *)malloc(xmax*sizeof(COLR));
277 if (scanln == NULL) {
278 perror(progname);
279 exit(1);
280 }
281 for (i = 0; i < 24; i++) {
282 setcolor(inpRGB[i], 0., 0., 0.);
283 ccount[i] = 0;
284 }
285 for (y = ymax-1; y >= 0; y--) {
286 if (freadcolrs(scanln, xmax, stdin) < 0) {
287 fprintf(stderr, "%s: error reading input picture\n",
288 progname);
289 exit(1);
290 }
291 for (x = 0; x < xmax; x++) {
292 i = chartndx(x, y, CENTCVG);
293 if (i >= 0) {
294 colr_color(pval, scanln[x]);
295 addcolor(inpRGB[i], pval);
296 ccount[i]++;
297 }
298 }
299 }
300 for (i = 0; i < 24; i++) { /* compute averages */
301 if (ccount[i] == 0)
302 continue;
303 d = 1./ccount[i];
304 scalecolor(inpRGB[i], d);
305 inpflags |= 1L<<i;
306 }
307 free((char *)scanln);
308 }
309
310
311 getcolors() /* get xyY colors from standard input */
312 {
313 int gotwhite = 0;
314 COLOR whiteclr;
315 int n;
316 float xyYin[3];
317
318 while (fgetval(stdin, 'i', &n) == 1) { /* read colors */
319 if (n < 0 | n > 24 ||
320 fgetval(stdin, 'f', &xyYin[0]) != 1 ||
321 fgetval(stdin, 'f', &xyYin[1]) != 1 ||
322 fgetval(stdin, 'f', &xyYin[2]) != 1 ||
323 xyYin[0] < 0. | xyYin[0] > 1. |
324 xyYin[1] < 0. | xyYin[1] > 1.) {
325 fprintf(stderr, "%s: bad color input data\n",
326 progname);
327 exit(1);
328 }
329 if (n == 0) { /* calibration white */
330 xyY2RGB(whiteclr, xyYin);
331 gotwhite++;
332 } else { /* standard color */
333 n--;
334 xyY2RGB(inpRGB[n], xyYin);
335 inpflags |= 1L<<n;
336 }
337 }
338 /* normalize colors */
339 if (!gotwhite) {
340 if (!(inpflags & 1L<<White)) {
341 fprintf(stderr, "%s: missing input for White\n",
342 progname);
343 exit(1);
344 }
345 setcolor(whiteclr,
346 colval(inpRGB[White],RED)/colval(mbRGB[White],RED),
347 colval(inpRGB[White],GRN)/colval(mbRGB[White],GRN),
348 colval(inpRGB[White],BLU)/colval(mbRGB[White],BLU));
349 }
350 for (n = 0; n < 24; n++)
351 if (inpflags & 1L<<n)
352 setcolor(inpRGB[n],
353 colval(inpRGB[n],RED)/colval(whiteclr,RED),
354 colval(inpRGB[n],GRN)/colval(whiteclr,GRN),
355 colval(inpRGB[n],BLU)/colval(whiteclr,BLU));
356 }
357
358
359 bresp(y, x) /* piecewise linear interpolation of primaries */
360 COLOR y, x;
361 {
362 double cv[3];
363 register int i, n;
364
365 for (i = 0; i < 3; i++) {
366 for (n = 0; n < NMBNEU-2; n++)
367 if (colval(x,i) < colval(bramp[n+1][0],i))
368 break;
369 cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
370 colval(bramp[n][1],i) +
371 (colval(x,i) - colval(bramp[n][0],i)) *
372 colval(bramp[n+1][1],i)) /
373 (colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
374 if (cv[i] < 0.) cv[i] = 0.;
375 }
376 setcolor(y, cv[0], cv[1], cv[2]);
377 }
378
379
380 compute() /* compute color mapping */
381 {
382 COLOR clrin[NMBMOD], clrout[NMBMOD];
383 register int i, n;
384 /* did we get what we need? */
385 if ((inpflags & REQFLGS) != REQFLGS) {
386 fprintf(stderr, "%s: missing required input colors\n",
387 progname);
388 exit(1);
389 }
390 /* compute piecewise luminance curve */
391 for (i = 0; i < NMBNEU; i++) {
392 copycolor(bramp[i][0], inpRGB[mbneu[i]]);
393 copycolor(bramp[i][1], mbRGB[mbneu[i]]);
394 }
395 /* compute color matrix */
396 for (n = 0, i = 0; i < NMBMOD; i++)
397 if (inpflags & 1L<<mbmod[i]) {
398 bresp(clrin[n], inpRGB[mbmod[i]]);
399 copycolor(clrout[n], mbRGB[mbmod[i]]);
400 n++;
401 }
402 compsoln(clrin, clrout, n);
403 }
404
405
406 putmapping() /* put out color mapping for pcomb -f */
407 {
408 static char cchar[3] = {'r', 'g', 'b'};
409 register int i, j;
410 /* print brightness mapping */
411 for (j = 0; j < 3; j++) {
412 printf("%cxa(i) : select(i", cchar[j]);
413 for (i = 0; i < NMBNEU; i++)
414 printf(",%g", colval(bramp[i][0],j));
415 printf(");\n");
416 printf("%cya(i) : select(i", cchar[j]);
417 for (i = 0; i < NMBNEU; i++)
418 printf(",%g", colval(bramp[i][1],j));
419 printf(");\n");
420 printf("%c = %ci(1);\n", cchar[j], cchar[j]);
421 printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
422 cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
423 cchar[j], cchar[j]);
424 printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
425 printf("%cn = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
426 cchar[j], cchar[j], cchar[j],
427 cchar[j], cchar[j], cchar[j]);
428 printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
429 cchar[j], cchar[j], cchar[j],
430 cchar[j], cchar[j]);
431 printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
432 cchar[j], cchar[j], cchar[j], cchar[j]);
433 }
434 /* print color mapping */
435 printf("ro = %g*rn + %g*gn + %g*bn ;\n",
436 solmat[0][0], solmat[0][1], solmat[0][2]);
437 printf("go = %g*rn + %g*gn + %g*bn ;\n",
438 solmat[1][0], solmat[1][1], solmat[1][2]);
439 printf("bo = %g*rn + %g*gn + %g*bn ;\n",
440 solmat[2][0], solmat[2][1], solmat[2][2]);
441 }
442
443
444 compsoln(cin, cout, n) /* solve 3xN system using least-squares */
445 COLOR cin[], cout[];
446 int n;
447 {
448 extern double mx3d_adjoint(), fabs();
449 double mat[3][3], invmat[3][3];
450 double det;
451 double colv[3], rowv[3];
452 register int i, j, k;
453
454 if (n < 3 | n > NMBMOD) {
455 fprintf(stderr, "%s: bad number of colors to match\n", progname);
456 exit(1);
457 }
458 if (n == 3)
459 for (i = 0; i < 3; i++)
460 for (j = 0; j < 3; j++)
461 mat[i][j] = colval(cin[j],i);
462 else { /* compute A^t A */
463 for (i = 0; i < 3; i++)
464 for (j = i; j < 3; j++) {
465 mat[i][j] = 0.;
466 for (k = 0; k < n; k++)
467 mat[i][j] += colval(cin[k],i) *
468 colval(cin[k],j);
469 }
470 for (i = 1; i < 3; i++) /* using symmetry */
471 for (j = 0; j < i; j++)
472 mat[i][j] = mat[j][i];
473 }
474 det = mx3d_adjoint(mat, invmat);
475 if (fabs(det) < 1e-4) {
476 fprintf(stderr, "%s: cannot compute color mapping\n",
477 progname);
478 solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
479 solmat[0][1] = solmat[0][2] = solmat[1][0] =
480 solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
481 return;
482 }
483 for (i = 0; i < 3; i++)
484 for (j = 0; j < 3; j++)
485 invmat[i][j] /= det;
486 for (i = 0; i < 3; i++) {
487 if (n == 3)
488 for (j = 0; j < 3; j++)
489 colv[j] = colval(cout[j],i);
490 else
491 for (j = 0; j < 3; j++) {
492 colv[j] = 0.;
493 for (k = 0; k < n; k++)
494 colv[j] += colval(cout[k],i) *
495 colval(cin[k],j);
496 }
497 mx3d_transform(colv, invmat, rowv);
498 for (j = 0; j < 3; j++)
499 solmat[i][j] = rowv[j];
500 }
501 }
502
503
504 cvtcolor(cout, cin) /* convert color according to our mapping */
505 COLOR cout, cin;
506 {
507 double r, g, b;
508
509 bresp(cout, cin);
510 r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1]
511 + colval(cout,2)*solmat[0][2];
512 if (r < 0) r = 0;
513 g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1]
514 + colval(cout,2)*solmat[1][2];
515 if (g < 0) g = 0;
516 b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
517 + colval(cout,2)*solmat[2][2];
518 if (b < 0) b = 0;
519 setcolor(cout, r, g, b);
520 }
521
522
523 xyY2RGB(rgbout, xyYin) /* convert xyY to RGB */
524 COLOR rgbout;
525 register float xyYin[3];
526 {
527 COLOR ctmp;
528 double d;
529
530 d = xyYin[2] / xyYin[1];
531 ctmp[0] = xyYin[0] * d;
532 ctmp[1] = xyYin[2];
533 ctmp[2] = (1. - xyYin[0] - xyYin[1]) * d;
534 cie_rgb(rgbout, ctmp);
535 }
536
537
538 picdebug() /* put out debugging picture */
539 {
540 COLOR *scan;
541 int y;
542 register int x, i;
543
544 if (fseek(stdin, 0L, 0) == EOF) {
545 fprintf(stderr, "%s: cannot seek on input picture\n", progname);
546 exit(1);
547 }
548 getheader(stdin, NULL, NULL); /* skip input header */
549 fgetresolu(&xmax, &ymax, stdin);
550 /* allocate scanline */
551 scan = (COLOR *)malloc(xmax*sizeof(COLOR));
552 if (scan == NULL) {
553 perror(progname);
554 exit(1);
555 }
556 /* finish debug header */
557 fputformat(COLRFMT, debugfp);
558 putc('\n', debugfp);
559 fprtresolu(xmax, ymax, debugfp);
560 /* write debug picture */
561 for (y = ymax-1; y >= 0; y--) {
562 if (freadscan(scan, xmax, stdin) < 0) {
563 fprintf(stderr, "%s: error rereading input picture\n",
564 progname);
565 exit(1);
566 }
567 for (x = 0; x < xmax; x++) {
568 i = chartndx(x, y, CENTCVG);
569 if (i < 0)
570 cvtcolor(scan[x], scan[x]);
571 else
572 copycolor(scan[x], mbRGB[i]);
573 }
574 if (fwritescan(scan, xmax, debugfp) < 0) {
575 fprintf(stderr, "%s: error writing debugging picture\n",
576 progname);
577 exit(1);
578 }
579 }
580 /* clean up */
581 fclose(debugfp);
582 free((char *)scan);
583 }
584
585
586 clrdebug() /* put out debug picture from color input */
587 {
588 static COLR blkclr = BLKCOLR;
589 COLR mbclr[24], cvclr[24];
590 COLR *scan;
591 COLOR ctmp;
592 int y;
593 register int i, x;
594 /* convert colors */
595 for (i = 0; i < 24; i++) {
596 setcolr(mbclr[i], colval(mbRGB[i],RED),
597 colval(mbRGB[i],GRN), colval(mbRGB[i],BLU));
598 if (inpflags & 1L<<i) {
599 cvtcolor(ctmp, inpRGB[i]);
600 setcolr(cvclr[i], colval(ctmp,RED),
601 colval(ctmp,GRN), colval(ctmp,BLU));
602 }
603 }
604 /* allocate scanline */
605 scan = (COLR *)malloc(xmax*sizeof(COLR));
606 if (scan == NULL) {
607 perror(progname);
608 exit(1);
609 }
610 /* finish debug header */
611 fputformat(COLRFMT, debugfp);
612 putc('\n', debugfp);
613 fprtresolu(xmax, ymax, debugfp);
614 /* write debug picture */
615 for (y = ymax-1; y >= 0; y--) {
616 for (x = 0; x < xmax; x++)
617 if ((i = chartndx(x, y, CENTCVG)) >= 0)
618 copycolr(scan[x], mbclr[i]);
619 else if ((i = chartndx(x, y, FULLCVG)) >= 0 &&
620 inpflags & 1L<<i)
621 copycolr(scan[x], cvclr[i]);
622 else
623 copycolr(scan[x], blkclr);
624 if (fwritecolrs(scan, xmax, debugfp) < 0) {
625 fprintf(stderr, "%s: error writing debugging picture\n",
626 progname);
627 exit(1);
628 }
629 }
630 /* clean up */
631 fclose(debugfp);
632 free((char *)scan);
633 }