ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/macbethcal.c
Revision: 2.3
Committed: Wed Oct 11 18:52:00 1995 UTC (28 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +63 -14 lines
Log Message:
non-working version using buggy svdcmp() routine

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     /* MacBeth colors (CIE 1931, absolute white) */
22     /* computed from spectral measurements */
23     float mbxyY[24][3] = {
24     {0.462, 0.3769, 0.0932961}, /* DarkSkin */
25     {0.4108, 0.3542, 0.410348}, /* LightSkin */
26     {0.2626, 0.267, 0.181554}, /* BlueSky */
27     {0.36, 0.4689, 0.108447}, /* Foliage */
28     {0.2977, 0.2602, 0.248407}, /* BlueFlower */
29     {0.2719, 0.3485, 0.401156}, /* BluishGreen */
30     {0.52, 0.4197, 0.357899}, /* Orange */
31     {0.229, 0.1866, 0.103911}, /* PurplishBlue */
32     {0.4909, 0.3262, 0.242615}, /* ModerateRed */
33     {0.3361, 0.2249, 0.0600102}, /* Purple */
34     {0.3855, 0.4874, 0.42963}, /* YellowGreen */
35     {0.4853, 0.4457, 0.476343}, /* OrangeYellow */
36     {0.2026, 0.1369, 0.0529249}, /* Blue */
37     {0.3007, 0.4822, 0.221226}, /* Green */
38     {0.5805, 0.3238, 0.162167}, /* Red */
39     {0.4617, 0.472, 0.64909}, /* Yellow */
40     {0.4178, 0.2625, 0.233662}, /* Magenta */
41     {0.2038, 0.2508, 0.167275}, /* Cyan */
42     {0.3358, 0.337, 0.916877}, /* White */
43     {0.3338, 0.3348, 0.604678}, /* Neutral.8 */
44     {0.3333, 0.3349, 0.364566}, /* Neutral.65 */
45     {0.3353, 0.3359, 0.200238}, /* Neutral.5 */
46     {0.3363, 0.336, 0.0878721}, /* Neutral.35 */
47     {0.3346, 0.3349, 0.0308383} /* Black */
48     };
49    
50     COLOR mbRGB[24]; /* MacBeth RGB values */
51    
52     #define NMBNEU 6 /* Number of MacBeth neutral colors */
53     short mbneu[NMBNEU] = {23,22,21,20,19,18};
54 greg 2.3 #define NMBMOD 16 /* Number of MacBeth unsaturated colors */
55     short mbmod[NMBMOD] = {0,1,2,3,4,5,6,7,8,9,10,11,19,20,21,22};
56 greg 2.1 #define NMBSAT 6 /* Number of MacBeth saturated colors */
57     short mbsat[NMBSAT] = {14,12,13,15,16,17};
58    
59     int xmax, ymax; /* input image dimensions */
60     int bounds[4][2]; /* image coordinates of chart corners */
61     double imgxfm[3][3]; /* coordinate transformation matrix */
62    
63     COLOR picRGB[24]; /* picture colors */
64    
65 greg 2.2 COLOR bramp[NMBNEU][2]; /* brightness ramp (per primary) */
66 greg 2.1 double solmat[3][3]; /* color mapping matrix */
67    
68     FILE *debugfp = NULL;
69     char *progname;
70    
71     extern char *malloc();
72    
73    
74     main(argc, argv)
75     int argc;
76     char **argv;
77     {
78     int i;
79    
80     progname = argv[0];
81     if (argc > 2 && !strcmp(argv[1], "-d")) { /* debug output */
82     if ((debugfp = fopen(argv[2], "w")) == NULL) {
83     perror(argv[2]);
84     exit(1);
85     }
86     #ifdef MSDOS
87     setmode(fileno(debugfp), O_BINARY);
88     #endif
89     newheader("RADIANCE", debugfp);
90     printargs(argc, argv, debugfp);
91     argv += 2;
92     argc -= 2;
93     }
94     if (argc != 3 && argc != 11)
95     goto userr;
96     if (strcmp(argv[1], "-") && freopen(argv[1], "r", stdin) == NULL) {
97     perror(argv[1]);
98     exit(1);
99     }
100     if (strcmp(argv[2], "-") && freopen(argv[2], "w", stdout) == NULL) {
101     perror(argv[2]);
102     exit(1);
103     }
104     #ifdef MSDOS
105     setmode(fileno(stdin), O_BINARY);
106     #endif
107     if (checkheader(stdin, COLRFMT, NULL) < 0 ||
108     fgetresolu(&xmax, &ymax, stdin) < 0) {
109     fprintf(stderr, "%s: bad input picture\n", progname);
110     exit(1);
111     }
112     /* get chart boundaries */
113     if (argc == 11) {
114     for (i = 0; i < 4; i++) {
115     if (!isint(argv[2*i+3]) | !isint(argv[2*i+4]))
116     goto userr;
117     bounds[i][0] = atoi(argv[2*i+3]);
118     bounds[i][1] = atoi(argv[2*i+4]);
119     }
120     } else {
121     bounds[0][0] = bounds[2][0] = .029*xmax + .5;
122     bounds[0][1] = bounds[1][1] = .956*ymax + .5;
123     bounds[1][0] = bounds[3][0] = .971*xmax + .5;
124     bounds[2][1] = bounds[3][1] = .056*ymax + .5;
125     }
126     init(); /* initialize */
127     getcolors(); /* get picture colors */
128 greg 2.2 compute(); /* compute color mapping */
129 greg 2.1 putmapping(); /* put out color mapping */
130     putdebug(); /* put out debug picture */
131     exit(0);
132     userr:
133     fprintf(stderr, "Usage: %s [-d dbg.pic] input.pic output.cal [xul yul xur yur xll yll xlr ylr]\n",
134     progname);
135     exit(1);
136     }
137    
138    
139     init() /* initialize */
140     {
141     double quad[4][2];
142     /* make coordinate transformation */
143     quad[0][0] = bounds[0][0];
144     quad[0][1] = bounds[0][1];
145     quad[1][0] = bounds[1][0];
146     quad[1][1] = bounds[1][1];
147     quad[2][0] = bounds[3][0];
148     quad[2][1] = bounds[3][1];
149     quad[3][0] = bounds[2][0];
150     quad[3][1] = bounds[2][1];
151    
152     if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
153     fprintf(stderr, "%s: bad chart boundaries\n", progname);
154     exit(1);
155     }
156     }
157    
158    
159     int
160     chartndx(x, y) /* find color number for position */
161     int x, y;
162     {
163     double ipos[3], cpos[3];
164     int ix, iy;
165     double fx, fy;
166    
167     ipos[0] = x;
168     ipos[1] = y;
169     ipos[2] = 1;
170     mx3d_transform(ipos, imgxfm, cpos);
171     cpos[0] /= cpos[2];
172     cpos[1] /= cpos[2];
173     if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
174     return(-1);
175     ix = cpos[0];
176     iy = cpos[1];
177     fx = cpos[0] - ix;
178     fy = cpos[1] - iy;
179     if (fx < .35 || fx >= .65 || fy < .35 || fy >= .65)
180     return(-1);
181     return(iy*6 + ix);
182     }
183    
184    
185     getcolors() /* load in picture colors */
186     {
187     COLR *scanln;
188     COLOR pval;
189     int ccount[24];
190     double d;
191     int y;
192     register int x, i;
193    
194     scanln = (COLR *)malloc(xmax*sizeof(COLR));
195     if (scanln == NULL) {
196     perror(progname);
197     exit(1);
198     }
199     for (i = 0; i < 24; i++) {
200     setcolor(picRGB[i], 0., 0., 0.);
201     ccount[i] = 0;
202     }
203     for (y = ymax-1; y >= 0; y--) {
204     if (freadcolrs(scanln, xmax, stdin) < 0) {
205     fprintf(stderr, "%s: error reading input picture\n",
206     progname);
207     exit(1);
208     }
209     for (x = 0; x < xmax; x++) {
210     i = chartndx(x, y);
211     if (i >= 0) {
212     colr_color(pval, scanln[x]);
213     addcolor(picRGB[i], pval);
214     ccount[i]++;
215     }
216     }
217     }
218     for (i = 0; i < 24; i++) {
219     if (ccount[i] == 0) {
220     fprintf(stderr, "%s: bad chart boundaries\n",
221     progname);
222     exit(1);
223     }
224     d = 1.0/ccount[i];
225     scalecolor(picRGB[i], d);
226     }
227     free((char *)scanln);
228     }
229    
230    
231 greg 2.2 bresp(y, x) /* piecewise linear interpolation of primaries */
232     COLOR y, x;
233 greg 2.1 {
234 greg 2.3 double cv[3];
235 greg 2.2 register int i, n;
236 greg 2.1
237 greg 2.2 for (i = 0; i < 3; i++) {
238     n = NMBNEU;
239     while (n > 0 && colval(x,i) < colval(bramp[--n][0],i))
240     ;
241 greg 2.3 cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
242 greg 2.2 colval(bramp[n][1],i) +
243     (colval(x,i) - colval(bramp[n][0],i)) *
244     colval(bramp[n+1][1],i)) /
245     (colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
246 greg 2.3 if (cv[i] < 0.) cv[i] = 0.;
247 greg 2.2 }
248 greg 2.3 setcolor(y, cv[0], cv[1], cv[2]);
249 greg 2.1 }
250    
251    
252 greg 2.2 compute() /* compute color mapping */
253 greg 2.1 {
254 greg 2.2 COLOR clrin[NMBMOD], clrout[NMBMOD];
255     COLOR ctmp;
256     double d;
257     register int i;
258     /* map MacBeth colors to RGB space */
259     for (i = 0; i < 24; i++) {
260     d = mbxyY[i][2] / mbxyY[i][1];
261     ctmp[0] = mbxyY[i][0] * d;
262     ctmp[1] = mbxyY[i][2];
263     ctmp[2] = (1. - mbxyY[i][0] - mbxyY[i][1]) * d;
264     cie_rgb(mbRGB[i], ctmp);
265     }
266 greg 2.1 /* compute piecewise luminance curve */
267     for (i = 0; i < NMBNEU; i++) {
268 greg 2.2 copycolor(bramp[i][0], picRGB[mbneu[i]]);
269     copycolor(bramp[i][1], mbRGB[mbneu[i]]);
270 greg 2.1 }
271     /* compute color matrix */
272 greg 2.2 for (i = 0; i < NMBMOD; i++) {
273     bresp(clrin[i], picRGB[mbmod[i]]);
274     copycolor(clrout[i], mbRGB[mbmod[i]]);
275     }
276 greg 2.1 compsoln(clrin, clrout, NMBMOD);
277 greg 2.2 }
278    
279    
280     putmapping() /* put out color mapping for pcomb -f */
281     {
282     static char cchar[3] = {'r', 'g', 'b'};
283     register int i, j;
284 greg 2.1 /* print brightness mapping */
285 greg 2.2 for (j = 0; j < 3; j++) {
286     printf("%cxa(i) : select(i", cchar[j]);
287     for (i = 0; i < NMBNEU; i++)
288     printf(",%g", colval(bramp[i][0],j));
289     printf(");\n");
290     printf("%cya(i) : select(i", cchar[j]);
291     for (i = 0; i < NMBNEU; i++)
292     printf(",%g", colval(bramp[i][1],j));
293     printf(");\n");
294     printf("%c = %ci(1);\n", cchar[j], cchar[j]);
295     printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
296     cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
297     cchar[j], cchar[j]);
298     printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
299     printf("%cn = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
300     cchar[j], cchar[j], cchar[j],
301     cchar[j], cchar[j], cchar[j]);
302     printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
303     cchar[j], cchar[j], cchar[j],
304     cchar[j], cchar[j]);
305     printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
306     cchar[j], cchar[j], cchar[j], cchar[j]);
307     }
308 greg 2.1 /* print color mapping */
309 greg 2.2 printf("ro = %g*rn + %g*gn + %g*bn ;\n",
310 greg 2.3 solmat[0][0], solmat[0][1], solmat[0][2]);
311 greg 2.2 printf("go = %g*rn + %g*gn + %g*bn ;\n",
312 greg 2.3 solmat[1][0], solmat[1][1], solmat[1][2]);
313 greg 2.2 printf("bo = %g*rn + %g*gn + %g*bn ;\n",
314 greg 2.3 solmat[2][0], solmat[2][1], solmat[2][2]);
315 greg 2.1 }
316    
317    
318 greg 2.3 #if NMBMOD == 3
319 greg 2.1 compsoln(cin, cout, n) /* solve 3x3 system */
320 greg 2.2 COLOR cin[], cout[];
321 greg 2.1 int n;
322     {
323     extern double mx3d_adjoint(), fabs();
324     double mat[3][3], invmat[3][3];
325     double det;
326     double colv[3], rowv[3];
327     register int i, j;
328    
329     if (n != 3) {
330     fprintf(stderr, "%s: inconsistent code!\n", progname);
331     exit(1);
332     }
333     for (i = 0; i < 3; i++)
334     for (j = 0; j < 3; j++)
335 greg 2.2 mat[i][j] = colval(cin[j],i);
336 greg 2.1 det = mx3d_adjoint(mat, invmat);
337     if (fabs(det) < 1e-4) {
338     fprintf(stderr, "%s: cannot compute color mapping\n",
339     progname);
340     solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
341     solmat[0][1] = solmat[0][2] = solmat[1][0] =
342     solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
343     return;
344     }
345     for (i = 0; i < 3; i++)
346     for (j = 0; j < 3; j++)
347     invmat[i][j] /= det;
348     for (i = 0; i < 3; i++) {
349     for (j = 0; j < 3; j++)
350 greg 2.3 colv[j] = colval(cout[j],i);
351     mx3d_transform(colv, invmat, rowv);
352 greg 2.1 for (j = 0; j < 3; j++)
353 greg 2.3 solmat[i][j] = rowv[j];
354 greg 2.1 }
355     }
356 greg 2.3 #else
357     compsoln(cin, cout, n) /* solve 3xN system (N > 3) */
358     COLOR cin[], cout[];
359     int n;
360     {
361     double *au[NMBMOD], *v[3], vv[3][3], auv[NMBMOD][3], w[3];
362     double b[NMBMOD];
363     register int i, j;
364 greg 2.1
365 greg 2.3 if (n > NMBMOD) {
366     fprintf(stderr, "%s: inconsistent code!\n", progname);
367     exit(1);
368     }
369     for (i = 0; i < n; i++) /* assign rectangular matrix A */
370     for (j = 0; j < 3; j++)
371     auv[i][j] = colval(cin[i],j);
372     /* svdcmp indexing requires pointer offsets */
373     for (j = 0; j < 3; j++)
374     v[j] = vv[j] - 1;
375     for (i = 0; i < n; i++)
376     au[i] = auv[i] - 1;
377     /* compute singular value decomposition */
378     fprintf(stderr, "A:\n");
379     for (i = 1; i <= n; i++)
380     fprintf(stderr, "%g %g %g\n", (au-1)[i][1], (au-1)[i][2], (au-1)[i][3]);
381     svdcmp(au-1, n, 3, w-1, v-1);
382     fprintf(stderr, "U:\n");
383     for (i = 0; i < n; i++)
384     fprintf(stderr, "%g %g %g\n", auv[i][0], auv[i][1], auv[i][2]);
385     fprintf(stderr, "V:\n");
386     for (i = 0; i < 3; i++)
387     fprintf(stderr, "%g %g %g\n", vv[i][0], vv[i][1], vv[i][2]);
388     fprintf(stderr, "W: %g %g %g\n", w[0], w[1], w[2]);
389     /* zero out small weights */
390     for (j = 0; j < 3; j++)
391     if (w[j] < 1e-4)
392     w[j] = 0.;
393     /* back substitution for each row vector */
394     for (j = 0; j < 3; j++) {
395     for (i = 0; i < n; i++)
396     b[i] = colval(cout[i],j);
397     svbksb(au-1, w-1, v-1, n, 3, b-1, solmat[j]-1);
398     }
399     }
400     #endif
401 greg 2.1
402 greg 2.3
403 greg 2.1 cvtcolor(cout, cin) /* convert color according to our mapping */
404     COLOR cout, cin;
405     {
406     double r, g, b;
407    
408 greg 2.2 bresp(cout, cin);
409 greg 2.3 r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1]
410     + colval(cout,2)*solmat[0][2];
411 greg 2.2 if (r < 0) r = 0;
412 greg 2.3 g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1]
413     + colval(cout,2)*solmat[1][2];
414 greg 2.2 if (g < 0) g = 0;
415 greg 2.3 b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
416 greg 2.2 + colval(cout,2)*solmat[2][2];
417     if (b < 0) b = 0;
418     setcolor(cout, r, g, b);
419 greg 2.1 }
420    
421    
422     putdebug() /* put out debugging picture */
423     {
424     COLOR *scan;
425     int y;
426     register int x, i;
427    
428     if (debugfp == NULL)
429     return;
430     if (fseek(stdin, 0L, 0) == EOF) {
431     fprintf(stderr, "%s: cannot seek on input picture\n", progname);
432     exit(1);
433     }
434     getheader(stdin, NULL, NULL); /* skip input header */
435     fgetresolu(&xmax, &ymax, stdin);
436     /* allocate scanline */
437     scan = (COLOR *)malloc(xmax*sizeof(COLOR));
438     if (scan == NULL) {
439     perror(progname);
440     exit(1);
441     }
442     /* finish debug header */
443     putc('\n', debugfp);
444     fprtresolu(xmax, ymax, debugfp);
445     for (y = ymax-1; y >= 0; y--) {
446     if (freadscan(scan, xmax, stdin) < 0) {
447     fprintf(stderr, "%s: error rereading input picture\n",
448     progname);
449     exit(1);
450     }
451     for (x = 0; x < xmax; x++) {
452     i = chartndx(x, y);
453     if (i < 0)
454     cvtcolor(scan[x], scan[x]);
455     else
456     copycolor(scan[x], mbRGB[i]);
457     }
458     if (fwritescan(scan, xmax, debugfp) < 0) {
459     fprintf(stderr, "%s: error writing debugging picture\n",
460     progname);
461     exit(1);
462     }
463     }
464     free((char *)scan);
465     }