ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/macbethcal.c
Revision: 2.11
Committed: Thu Jan 30 17:07:42 1997 UTC (27 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.10: +41 -22 lines
Log Message:
changed ordering so color trans. before non-lin. corr. for output

File Contents

# User Rev Content
1 greg 2.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 greg 2.4 /* 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 greg 2.1 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 greg 2.4 short mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};
79 greg 2.1
80 greg 2.8 #define NEUFLGS (1L<<White|1L<<Neutral8|1L<<Neutral65| \
81     1L<<Neutral5|1L<<Neutral35|1L<<Black)
82 greg 2.4
83 greg 2.8 #define SATFLGS (1L<<Red|1L<<Green|1L<<Blue|1L<<Magenta|1L<<Yellow| \
84     1L<<Cyan|1L<<Orange|1L<<Purple|1L<<PurplishBlue| \
85     1L<<YellowGreen|1<<OrangeYellow|1L<<BlueFlower)
86 greg 2.4
87 greg 2.8 #define UNSFLGS (1L<<DarkSkin|1L<<LightSkin|1L<<BlueSky|1L<<Foliage| \
88     1L<<BluishGreen|1L<<ModerateRed)
89 greg 2.7
90 greg 2.8 #define REQFLGS NEUFLGS /* need these colors */
91     #define MODFLGS (NEUFLGS|UNSFLGS) /* should be in gamut */
92    
93 greg 2.9 #define RG_BORD 0 /* patch border */
94     #define RG_CENT 01 /* central region of patch */
95     #define RG_ORIG 02 /* original color region */
96     #define RG_CORR 04 /* corrected color region */
97 greg 2.7
98 greg 2.11 int scanning = 1; /* scanned input (or recorded output)? */
99    
100 greg 2.1 int xmax, ymax; /* input image dimensions */
101     int bounds[4][2]; /* image coordinates of chart corners */
102     double imgxfm[3][3]; /* coordinate transformation matrix */
103    
104 greg 2.7 COLOR inpRGB[24]; /* measured or scanned input colors */
105     long inpflags = 0; /* flags of which colors were input */
106 greg 2.8 long gmtflags = 0; /* flags of out-of-gamut colors */
107 greg 2.1
108 greg 2.2 COLOR bramp[NMBNEU][2]; /* brightness ramp (per primary) */
109 greg 2.1 double solmat[3][3]; /* color mapping matrix */
110    
111 greg 2.4 FILE *debugfp = NULL; /* debug output picture */
112 greg 2.1 char *progname;
113    
114     extern char *malloc();
115    
116    
117     main(argc, argv)
118     int argc;
119     char **argv;
120     {
121     int i;
122    
123     progname = argv[0];
124 greg 2.7 for (i = 1; i < argc && argv[i][0] == '-'; i++)
125     switch (argv[i][1]) {
126     case 'd': /* debug output */
127     i++;
128     if (badarg(argc-i, argv+i, "s"))
129     goto userr;
130     if ((debugfp = fopen(argv[i], "w")) == NULL) {
131     perror(argv[i]);
132     exit(1);
133     }
134 greg 2.1 #ifdef MSDOS
135 greg 2.7 setmode(fileno(debugfp), O_BINARY);
136 greg 2.1 #endif
137 greg 2.7 newheader("RADIANCE", debugfp); /* start */
138     printargs(argc, argv, debugfp); /* header */
139     break;
140     case 'p': /* picture position */
141     if (badarg(argc-i-1, argv+i+1, "iiiiiiii"))
142     goto userr;
143     bounds[0][0] = atoi(argv[++i]);
144     bounds[0][1] = atoi(argv[++i]);
145     bounds[1][0] = atoi(argv[++i]);
146     bounds[1][1] = atoi(argv[++i]);
147     bounds[2][0] = atoi(argv[++i]);
148     bounds[2][1] = atoi(argv[++i]);
149     bounds[3][0] = atoi(argv[++i]);
150     bounds[3][1] = atoi(argv[++i]);
151 greg 2.11 scanning = 2;
152 greg 2.7 break;
153     case 'c': /* color input */
154 greg 2.11 scanning = 0;
155 greg 2.7 break;
156     default:
157     goto userr;
158     }
159     /* open files */
160     if (i < argc && freopen(argv[i], "r", stdin) == NULL) {
161 greg 2.1 perror(argv[1]);
162     exit(1);
163     }
164 greg 2.7 if (i+1 < argc && freopen(argv[i+1], "w", stdout) == NULL) {
165 greg 2.1 perror(argv[2]);
166     exit(1);
167     }
168 greg 2.11 if (scanning) { /* load input picture header */
169 greg 2.1 #ifdef MSDOS
170 greg 2.7 setmode(fileno(stdin), O_BINARY);
171 greg 2.1 #endif
172 greg 2.7 if (checkheader(stdin, COLRFMT, NULL) < 0 ||
173     fgetresolu(&xmax, &ymax, stdin) < 0) {
174     fprintf(stderr, "%s: bad input picture\n", progname);
175     exit(1);
176     }
177     } else { /* else set default xmax and ymax */
178     xmax = 512;
179     ymax = 2*512/3;
180 greg 2.1 }
181 greg 2.11 if (scanning != 2) { /* use default boundaries */
182 greg 2.1 bounds[0][0] = bounds[2][0] = .029*xmax + .5;
183     bounds[0][1] = bounds[1][1] = .956*ymax + .5;
184     bounds[1][0] = bounds[3][0] = .971*xmax + .5;
185     bounds[2][1] = bounds[3][1] = .056*ymax + .5;
186     }
187     init(); /* initialize */
188 greg 2.11 if (scanning) /* get picture colors */
189 greg 2.7 getpicture();
190     else
191     getcolors();
192 greg 2.2 compute(); /* compute color mapping */
193 greg 2.4 /* print comment */
194 greg 2.10 printf("{\n\tColor correction file computed by:\n\t\t");
195 greg 2.7 printargs(argc, argv, stdout);
196 greg 2.10 printf("\n\tUsage: pcomb -f %s uncorrected.pic > corrected.pic\n",
197     i+1 < argc ? argv[i+1] : "{this_file}");
198 greg 2.11 if (!scanning)
199     printf("\t Or: pcond [options] -f %s orig.pic > output.pic\n",
200     i+1 < argc ? argv[i+1] : "{this_file}");
201 greg 2.7 printf("}\n");
202 greg 2.1 putmapping(); /* put out color mapping */
203 greg 2.7 if (debugfp != NULL) /* put out debug picture */
204 greg 2.11 if (scanning)
205 greg 2.7 picdebug();
206     else
207     clrdebug();
208 greg 2.1 exit(0);
209     userr:
210 greg 2.7 fprintf(stderr,
211     "Usage: %s [-d dbg.pic][-p xul yul xur yur xll yll xlr ylr] input.pic [output.cal]\n",
212 greg 2.1 progname);
213 greg 2.7 fprintf(stderr, " or: %s [-d dbg.pic] -c [xyY.dat [output.cal]]\n",
214     progname);
215 greg 2.1 exit(1);
216     }
217    
218    
219     init() /* initialize */
220     {
221     double quad[4][2];
222 greg 2.7 register int i;
223 greg 2.1 /* make coordinate transformation */
224     quad[0][0] = bounds[0][0];
225     quad[0][1] = bounds[0][1];
226     quad[1][0] = bounds[1][0];
227     quad[1][1] = bounds[1][1];
228     quad[2][0] = bounds[3][0];
229     quad[2][1] = bounds[3][1];
230     quad[3][0] = bounds[2][0];
231     quad[3][1] = bounds[2][1];
232    
233     if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
234     fprintf(stderr, "%s: bad chart boundaries\n", progname);
235     exit(1);
236     }
237 greg 2.7 /* map MacBeth colors to RGB space */
238     for (i = 0; i < 24; i++)
239     xyY2RGB(mbRGB[i], mbxyY[i]);
240 greg 2.1 }
241    
242    
243     int
244 greg 2.9 chartndx(x, y, np) /* find color number for position */
245 greg 2.1 int x, y;
246 greg 2.9 int *np;
247 greg 2.1 {
248     double ipos[3], cpos[3];
249     int ix, iy;
250     double fx, fy;
251    
252     ipos[0] = x;
253     ipos[1] = y;
254     ipos[2] = 1;
255     mx3d_transform(ipos, imgxfm, cpos);
256     cpos[0] /= cpos[2];
257     cpos[1] /= cpos[2];
258     if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
259 greg 2.9 return(RG_BORD);
260 greg 2.1 ix = cpos[0];
261     iy = cpos[1];
262     fx = cpos[0] - ix;
263     fy = cpos[1] - iy;
264 greg 2.9 *np = iy*6 + ix;
265     if (fx >= 0.35 && fx < 0.65 && fy >= 0.35 && fy < 0.65)
266     return(RG_CENT);
267     if (fx < 0.05 || fx >= 0.95 || fy < 0.05 || fy >= 0.95)
268     return(RG_BORD);
269     if (fx >= 0.5) /* right side is corrected */
270     return(RG_CORR);
271     return(RG_ORIG); /* left side is original */
272 greg 2.1 }
273    
274    
275 greg 2.7 getpicture() /* load in picture colors */
276 greg 2.1 {
277     COLR *scanln;
278     COLOR pval;
279     int ccount[24];
280     double d;
281 greg 2.9 int y, i;
282     register int x;
283 greg 2.1
284     scanln = (COLR *)malloc(xmax*sizeof(COLR));
285     if (scanln == NULL) {
286     perror(progname);
287     exit(1);
288     }
289     for (i = 0; i < 24; i++) {
290 greg 2.7 setcolor(inpRGB[i], 0., 0., 0.);
291 greg 2.1 ccount[i] = 0;
292     }
293     for (y = ymax-1; y >= 0; y--) {
294     if (freadcolrs(scanln, xmax, stdin) < 0) {
295     fprintf(stderr, "%s: error reading input picture\n",
296     progname);
297     exit(1);
298     }
299 greg 2.9 for (x = 0; x < xmax; x++)
300     if (chartndx(x, y, &i) == RG_CENT) {
301 greg 2.1 colr_color(pval, scanln[x]);
302 greg 2.7 addcolor(inpRGB[i], pval);
303 greg 2.1 ccount[i]++;
304     }
305     }
306 greg 2.7 for (i = 0; i < 24; i++) { /* compute averages */
307     if (ccount[i] == 0)
308     continue;
309     d = 1./ccount[i];
310     scalecolor(inpRGB[i], d);
311     inpflags |= 1L<<i;
312     }
313     free((char *)scanln);
314     }
315    
316    
317     getcolors() /* get xyY colors from standard input */
318     {
319     int gotwhite = 0;
320     COLOR whiteclr;
321     int n;
322     float xyYin[3];
323    
324     while (fgetval(stdin, 'i', &n) == 1) { /* read colors */
325     if (n < 0 | n > 24 ||
326     fgetval(stdin, 'f', &xyYin[0]) != 1 ||
327     fgetval(stdin, 'f', &xyYin[1]) != 1 ||
328     fgetval(stdin, 'f', &xyYin[2]) != 1 ||
329     xyYin[0] < 0. | xyYin[0] > 1. |
330     xyYin[1] < 0. | xyYin[1] > 1.) {
331     fprintf(stderr, "%s: bad color input data\n",
332 greg 2.1 progname);
333     exit(1);
334     }
335 greg 2.7 if (n == 0) { /* calibration white */
336     xyY2RGB(whiteclr, xyYin);
337     gotwhite++;
338     } else { /* standard color */
339     n--;
340     xyY2RGB(inpRGB[n], xyYin);
341     inpflags |= 1L<<n;
342     }
343 greg 2.1 }
344 greg 2.7 /* normalize colors */
345     if (!gotwhite) {
346     if (!(inpflags & 1L<<White)) {
347     fprintf(stderr, "%s: missing input for White\n",
348     progname);
349     exit(1);
350     }
351     setcolor(whiteclr,
352     colval(inpRGB[White],RED)/colval(mbRGB[White],RED),
353     colval(inpRGB[White],GRN)/colval(mbRGB[White],GRN),
354     colval(inpRGB[White],BLU)/colval(mbRGB[White],BLU));
355     }
356     for (n = 0; n < 24; n++)
357     if (inpflags & 1L<<n)
358     setcolor(inpRGB[n],
359     colval(inpRGB[n],RED)/colval(whiteclr,RED),
360     colval(inpRGB[n],GRN)/colval(whiteclr,GRN),
361     colval(inpRGB[n],BLU)/colval(whiteclr,BLU));
362 greg 2.1 }
363    
364    
365 greg 2.2 bresp(y, x) /* piecewise linear interpolation of primaries */
366     COLOR y, x;
367 greg 2.1 {
368 greg 2.2 register int i, n;
369 greg 2.1
370 greg 2.2 for (i = 0; i < 3; i++) {
371 greg 2.5 for (n = 0; n < NMBNEU-2; n++)
372     if (colval(x,i) < colval(bramp[n+1][0],i))
373     break;
374 greg 2.8 colval(y,i) = ((colval(bramp[n+1][0],i) - colval(x,i)) *
375 greg 2.2 colval(bramp[n][1],i) +
376     (colval(x,i) - colval(bramp[n][0],i)) *
377     colval(bramp[n+1][1],i)) /
378     (colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
379     }
380 greg 2.1 }
381    
382    
383 greg 2.2 compute() /* compute color mapping */
384 greg 2.1 {
385 greg 2.8 COLOR clrin[24], clrout[24];
386     long cflags;
387     COLOR ctmp;
388     register int i, j, n;
389 greg 2.7 /* did we get what we need? */
390     if ((inpflags & REQFLGS) != REQFLGS) {
391     fprintf(stderr, "%s: missing required input colors\n",
392     progname);
393     exit(1);
394 greg 2.2 }
395 greg 2.1 /* compute piecewise luminance curve */
396     for (i = 0; i < NMBNEU; i++) {
397 greg 2.7 copycolor(bramp[i][0], inpRGB[mbneu[i]]);
398 greg 2.2 copycolor(bramp[i][1], mbRGB[mbneu[i]]);
399 greg 2.1 }
400 greg 2.8 /* compute color mapping */
401     do {
402     cflags = inpflags & ~gmtflags;
403     n = 0; /* compute transform matrix */
404     for (i = 0; i < 24; i++)
405     if (cflags & 1L<<i) {
406     bresp(clrin[n], inpRGB[i]);
407     copycolor(clrout[n], mbRGB[i]);
408     n++;
409     }
410     compsoln(clrin, clrout, n);
411     /* check out-of-gamut colors */
412     for (i = 0; i < 24; i++)
413     if (cflags & 1L<<i) {
414 greg 2.11 cvtcolor(ctmp, mbRGB[i]);
415 greg 2.8 for (j = 0; j < 3; j++)
416 greg 2.11 if (colval(ctmp,j) <= 1e-6 ||
417     colval(ctmp,j) >= 1.-1e-6) {
418 greg 2.8 gmtflags |= 1L<<i;
419     break;
420     }
421     }
422     } while (cflags & gmtflags);
423     if (gmtflags & MODFLGS)
424     fprintf(stderr,
425     "%s: warning - some moderate colors are out of gamut\n",
426     progname);
427 greg 2.2 }
428    
429    
430     putmapping() /* put out color mapping for pcomb -f */
431     {
432     static char cchar[3] = {'r', 'g', 'b'};
433     register int i, j;
434 greg 2.1 /* print brightness mapping */
435 greg 2.2 for (j = 0; j < 3; j++) {
436     printf("%cxa(i) : select(i", cchar[j]);
437     for (i = 0; i < NMBNEU; i++)
438     printf(",%g", colval(bramp[i][0],j));
439     printf(");\n");
440     printf("%cya(i) : select(i", cchar[j]);
441     for (i = 0; i < NMBNEU; i++)
442     printf(",%g", colval(bramp[i][1],j));
443     printf(");\n");
444     printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
445     cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
446     cchar[j], cchar[j]);
447     printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
448 greg 2.11 printf("%c%c = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
449     cchar[j], scanning?'n':'o', cchar[j],
450     cchar[j], cchar[j], cchar[j], cchar[j]);
451 greg 2.2 printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
452     cchar[j], cchar[j], cchar[j],
453     cchar[j], cchar[j]);
454     printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
455     cchar[j], cchar[j], cchar[j], cchar[j]);
456     }
457 greg 2.1 /* print color mapping */
458 greg 2.11 if (scanning) {
459     printf("r = ri(1); g = gi(1); b = bi(1);\n");
460     printf("ro = %g*rn + %g*gn + %g*bn ;\n",
461     solmat[0][0], solmat[0][1], solmat[0][2]);
462     printf("go = %g*rn + %g*gn + %g*bn ;\n",
463     solmat[1][0], solmat[1][1], solmat[1][2]);
464     printf("bo = %g*rn + %g*gn + %g*bn ;\n",
465     solmat[2][0], solmat[2][1], solmat[2][2]);
466     } else {
467     printf("r1 = ri(1); g1 = gi(1); b1 = bi(1);\n");
468     printf("r = %g*r1 + %g*g1 + %g*b1 ;\n",
469     solmat[0][0], solmat[0][1], solmat[0][2]);
470     printf("g = %g*r1 + %g*g1 + %g*b1 ;\n",
471     solmat[1][0], solmat[1][1], solmat[1][2]);
472     printf("b = %g*r1 + %g*g1 + %g*b1 ;\n",
473     solmat[2][0], solmat[2][1], solmat[2][2]);
474     }
475 greg 2.1 }
476    
477    
478 greg 2.4 compsoln(cin, cout, n) /* solve 3xN system using least-squares */
479 greg 2.2 COLOR cin[], cout[];
480 greg 2.1 int n;
481     {
482     extern double mx3d_adjoint(), fabs();
483     double mat[3][3], invmat[3][3];
484     double det;
485     double colv[3], rowv[3];
486 greg 2.4 register int i, j, k;
487 greg 2.1
488 greg 2.8 if (n < 3) {
489     fprintf(stderr, "%s: too few colors to match!\n", progname);
490 greg 2.1 exit(1);
491     }
492 greg 2.4 if (n == 3)
493     for (i = 0; i < 3; i++)
494     for (j = 0; j < 3; j++)
495     mat[i][j] = colval(cin[j],i);
496     else { /* compute A^t A */
497     for (i = 0; i < 3; i++)
498     for (j = i; j < 3; j++) {
499     mat[i][j] = 0.;
500     for (k = 0; k < n; k++)
501     mat[i][j] += colval(cin[k],i) *
502     colval(cin[k],j);
503     }
504     for (i = 1; i < 3; i++) /* using symmetry */
505     for (j = 0; j < i; j++)
506     mat[i][j] = mat[j][i];
507     }
508 greg 2.1 det = mx3d_adjoint(mat, invmat);
509     if (fabs(det) < 1e-4) {
510     fprintf(stderr, "%s: cannot compute color mapping\n",
511     progname);
512     solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
513     solmat[0][1] = solmat[0][2] = solmat[1][0] =
514     solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
515     return;
516     }
517     for (i = 0; i < 3; i++)
518     for (j = 0; j < 3; j++)
519     invmat[i][j] /= det;
520     for (i = 0; i < 3; i++) {
521 greg 2.4 if (n == 3)
522     for (j = 0; j < 3; j++)
523     colv[j] = colval(cout[j],i);
524     else
525     for (j = 0; j < 3; j++) {
526     colv[j] = 0.;
527     for (k = 0; k < n; k++)
528     colv[j] += colval(cout[k],i) *
529     colval(cin[k],j);
530     }
531 greg 2.3 mx3d_transform(colv, invmat, rowv);
532 greg 2.1 for (j = 0; j < 3; j++)
533 greg 2.3 solmat[i][j] = rowv[j];
534 greg 2.1 }
535     }
536    
537 greg 2.3
538 greg 2.1 cvtcolor(cout, cin) /* convert color according to our mapping */
539     COLOR cout, cin;
540     {
541 greg 2.8 COLOR ctmp;
542    
543 greg 2.11 if (scanning) {
544     bresp(ctmp, cin);
545     cresp(cout, ctmp);
546     } else {
547     cresp(ctmp, cin);
548     bresp(cout, ctmp);
549     }
550 greg 2.8 if (colval(cout,RED) < 0.)
551     colval(cout,RED) = 0.;
552     if (colval(cout,GRN) < 0.)
553     colval(cout,GRN) = 0.;
554     if (colval(cout,BLU) < 0.)
555     colval(cout,BLU) = 0.;
556     }
557    
558    
559     cresp(cout, cin) /* transform color according to matrix */
560     COLOR cout, cin;
561     {
562 greg 2.1 double r, g, b;
563    
564 greg 2.8 r = colval(cin,0)*solmat[0][0] + colval(cin,1)*solmat[0][1]
565     + colval(cin,2)*solmat[0][2];
566     g = colval(cin,0)*solmat[1][0] + colval(cin,1)*solmat[1][1]
567     + colval(cin,2)*solmat[1][2];
568     b = colval(cin,0)*solmat[2][0] + colval(cin,1)*solmat[2][1]
569     + colval(cin,2)*solmat[2][2];
570 greg 2.2 setcolor(cout, r, g, b);
571 greg 2.1 }
572    
573    
574 greg 2.7 xyY2RGB(rgbout, xyYin) /* convert xyY to RGB */
575     COLOR rgbout;
576     register float xyYin[3];
577 greg 2.1 {
578 greg 2.7 COLOR ctmp;
579     double d;
580    
581     d = xyYin[2] / xyYin[1];
582     ctmp[0] = xyYin[0] * d;
583     ctmp[1] = xyYin[2];
584     ctmp[2] = (1. - xyYin[0] - xyYin[1]) * d;
585     cie_rgb(rgbout, ctmp);
586     }
587    
588    
589     picdebug() /* put out debugging picture */
590     {
591 greg 2.8 static COLOR blkcol = BLKCOLOR;
592 greg 2.1 COLOR *scan;
593 greg 2.9 int y, i;
594     register int x, rg;
595 greg 2.1
596     if (fseek(stdin, 0L, 0) == EOF) {
597     fprintf(stderr, "%s: cannot seek on input picture\n", progname);
598     exit(1);
599     }
600     getheader(stdin, NULL, NULL); /* skip input header */
601     fgetresolu(&xmax, &ymax, stdin);
602     /* allocate scanline */
603     scan = (COLOR *)malloc(xmax*sizeof(COLOR));
604     if (scan == NULL) {
605     perror(progname);
606     exit(1);
607     }
608     /* finish debug header */
609 greg 2.5 fputformat(COLRFMT, debugfp);
610 greg 2.1 putc('\n', debugfp);
611     fprtresolu(xmax, ymax, debugfp);
612 greg 2.7 /* write debug picture */
613 greg 2.1 for (y = ymax-1; y >= 0; y--) {
614     if (freadscan(scan, xmax, stdin) < 0) {
615     fprintf(stderr, "%s: error rereading input picture\n",
616     progname);
617     exit(1);
618     }
619     for (x = 0; x < xmax; x++) {
620 greg 2.9 rg = chartndx(x, y, &i);
621     if (rg == RG_CENT) {
622     if (!(1L<<i & gmtflags) || (x+y)&07)
623     copycolor(scan[x], mbRGB[i]);
624     else
625     copycolor(scan[x], blkcol);
626     } else if (rg == RG_CORR)
627 greg 2.1 cvtcolor(scan[x], scan[x]);
628 greg 2.9 else if (rg != RG_ORIG)
629 greg 2.8 copycolor(scan[x], blkcol);
630 greg 2.1 }
631     if (fwritescan(scan, xmax, debugfp) < 0) {
632     fprintf(stderr, "%s: error writing debugging picture\n",
633     progname);
634     exit(1);
635     }
636     }
637 greg 2.7 /* clean up */
638     fclose(debugfp);
639     free((char *)scan);
640     }
641    
642    
643     clrdebug() /* put out debug picture from color input */
644     {
645     static COLR blkclr = BLKCOLR;
646 greg 2.9 COLR mbclr[24], cvclr[24], orclr[24];
647 greg 2.7 COLR *scan;
648     COLOR ctmp;
649 greg 2.9 int y, i;
650     register int x, rg;
651 greg 2.7 /* convert colors */
652     for (i = 0; i < 24; i++) {
653     setcolr(mbclr[i], colval(mbRGB[i],RED),
654     colval(mbRGB[i],GRN), colval(mbRGB[i],BLU));
655     if (inpflags & 1L<<i) {
656 greg 2.9 setcolr(orclr[i], colval(inpRGB[i],RED),
657     colval(inpRGB[i],GRN),
658     colval(inpRGB[i],BLU));
659 greg 2.7 cvtcolor(ctmp, inpRGB[i]);
660     setcolr(cvclr[i], colval(ctmp,RED),
661     colval(ctmp,GRN), colval(ctmp,BLU));
662     }
663     }
664     /* allocate scanline */
665     scan = (COLR *)malloc(xmax*sizeof(COLR));
666     if (scan == NULL) {
667     perror(progname);
668     exit(1);
669     }
670     /* finish debug header */
671     fputformat(COLRFMT, debugfp);
672     putc('\n', debugfp);
673     fprtresolu(xmax, ymax, debugfp);
674     /* write debug picture */
675     for (y = ymax-1; y >= 0; y--) {
676 greg 2.9 for (x = 0; x < xmax; x++) {
677     rg = chartndx(x, y, &i);
678     if (rg == RG_CENT) {
679 greg 2.8 if (!(1L<<i & gmtflags) || (x+y)&07)
680     copycolr(scan[x], mbclr[i]);
681     else
682     copycolr(scan[x], blkclr);
683 greg 2.9 } else if (rg == RG_BORD || !(1L<<i & inpflags))
684     copycolr(scan[x], blkclr);
685     else if (rg == RG_ORIG)
686     copycolr(scan[x], orclr[i]);
687     else /* rg == RG_CORR */
688 greg 2.7 copycolr(scan[x], cvclr[i]);
689 greg 2.9 }
690 greg 2.7 if (fwritecolrs(scan, xmax, debugfp) < 0) {
691     fprintf(stderr, "%s: error writing debugging picture\n",
692     progname);
693     exit(1);
694     }
695     }
696     /* clean up */
697     fclose(debugfp);
698 greg 2.1 free((char *)scan);
699     }