ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/px/macbethcal.c

Comparing ray/src/px/macbethcal.c (file contents):
Revision 2.1 by greg, Wed Oct 11 10:40:13 1995 UTC vs.
Revision 2.5 by greg, Tue Oct 17 10:24:21 1995 UTC

#	Line 18 | Line 18 | static char SCCSid[] = "$SunId$ LBL";
18		#include "resolu.h"
19		#include "pmap.h"
20
21	<	/* MacBeth colors (CIE 1931, absolute white) */
22	<	/* computed from spectral measurements */
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 */
#	Line 50 | Line 75 | float  mbxyY[24][3] = {
75		COLOR   mbRGB[24];              /* MacBeth RGB values */
76
77		#define NMBNEU          6       /* Number of MacBeth neutral colors */
78	<	short   mbneu[NMBNEU] = {23,22,21,20,19,18};
54	<	#define NMBMOD          3       /* Number of MacBeth moderate colors */
55	<	short   mbmod[NMBMOD] = {1,3,2};
56	<	#define NMBSAT          6       /* Number of MacBeth saturated colors */
57	<	short   mbsat[NMBSAT] = {14,12,13,15,16,17};
78	>	short   mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};
79
80	+	#define NMBMOD          17      /* Number of MacBeth unsaturated colors */
81	+	short   mbmod[NMBMOD] = {
82	+	DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen,
83	+	Orange,PurplishBlue,ModerateRed,Purple,YellowGreen,
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,OrangeYellow
91	+	};
92	+
93		int     xmax, ymax;             /* input image dimensions */
94		int     bounds[4][2];           /* image coordinates of chart corners */
95		double  imgxfm[3][3];           /* coordinate transformation matrix */
96
97		COLOR   picRGB[24];             /* picture colors */
98
99	<	double  bramp[NMBNEU][2];       /* brightness ramp */
99	>	COLOR   bramp[NMBNEU][2];       /* brightness ramp (per primary) */
100		double  solmat[3][3];           /* color mapping matrix */
101
102	<	FILE    *debugfp = NULL;
102	>	FILE    *debugfp = NULL;        /* debug output picture */
103		char    *progname;
104
105		extern char     *malloc();
#	Line 125 | Line 159 | char   **argv;
159		}
160		init();                         /* initialize */
161		getcolors();                    /* get picture colors */
162	+	compute();                      /* compute color mapping */
163	+	/* print comment */
164	+	printf("{ Color correction file computed by %s }\n", progname);
165	+	printf("{ from scanned MacBetch color chart %s }\n", argv[1]);
166		putmapping();                   /* put out color mapping */
167		putdebug();                     /* put out debug picture */
168		exit(0);
#	Line 138 | Line 176 | userr:
176		init()                          /* initialize */
177		{
178		double  quad[4][2];
141	–	COLOR   ctmp;
142	–	double  d;
143	–	int     i;
179		/* make coordinate transformation */
180		quad[0][0] = bounds[0][0];
181		quad[0][1] = bounds[0][1];
#	Line 155 | Line 190 | init()                         /* initialize */
190		fprintf(stderr, "%s: bad chart boundaries\n", progname);
191		exit(1);
192		}
158	–	/* map MacBeth colors to RGB space */
159	–	for (i = 0; i < 24; i++) {
160	–	d = mbxyY[i][2] / mbxyY[i][1];
161	–	ctmp[0] = mbxyY[i][0] * d;
162	–	ctmp[1] = mbxyY[i][2];
163	–	ctmp[2] = (1. - mbxyY[i][0] - mbxyY[i][1]) * d;
164	–	cie_rgb(mbRGB[i], ctmp);
165	–	}
193		}
194
195
#	Line 238 | Line 265 | getcolors()                            /* load in picture colors */
265		}
266
267
268	<	double
269	<	bresp(x)                /* piecewise linear interpolation of brightness */
243	<	double  x;
268	>	bresp(y, x)             /* piecewise linear interpolation of primaries */
269	>	COLOR   y, x;
270		{
271	<	register int    n = NMBNEU;
271	>	double  cv[3];
272	>	register int    i, n;
273
274	<	while (n > 0 && x < bramp[--n][0])
275	<	;
276	<	return( ((bramp[n+1][0] - x)*bramp[n][1] +
277	<	(x - bramp[n][0])*bramp[n+1][1]) /
278	<	(bramp[n+1][0] - bramp[n][0]) );
274	>	for (i = 0; i < 3; i++) {
275	>	for (n = 0; n < NMBNEU-2; n++)
276	>	if (colval(x,i) < colval(bramp[n+1][0],i))
277	>	break;
278	>	cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
279	>	colval(bramp[n][1],i) +
280	>	(colval(x,i) - colval(bramp[n][0],i)) *
281	>	colval(bramp[n+1][1],i)) /
282	>	(colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
283	>	if (cv[i] < 0.) cv[i] = 0.;
284	>	}
285	>	setcolor(y, cv[0], cv[1], cv[2]);
286		}
287
288
289	<	putmapping()                    /* compute and print mapping for pcomb -f */
289	>	compute()                       /* compute color mapping */
290		{
291	<	float   clrin[NMBMOD][3], clrout[NMBMOD][3];
292	<	register int    i, j;
291	>	COLOR   clrin[NMBMOD], clrout[NMBMOD];
292	>	COLOR   ctmp;
293	>	double  d;
294	>	register int    i;
295	>	/* map MacBeth colors to RGB space */
296	>	for (i = 0; i < 24; i++) {
297	>	d = mbxyY[i][2] / mbxyY[i][1];
298	>	ctmp[0] = mbxyY[i][0] * d;
299	>	ctmp[1] = mbxyY[i][2];
300	>	ctmp[2] = (1. - mbxyY[i][0] - mbxyY[i][1]) * d;
301	>	cie_rgb(mbRGB[i], ctmp);
302	>	}
303		/* compute piecewise luminance curve */
304		for (i = 0; i < NMBNEU; i++) {
305	<	bramp[i][0] = bright(picRGB[mbneu[i]]);
306	<	bramp[i][1] = bright(mbRGB[mbneu[i]]);
305	>	copycolor(bramp[i][0], picRGB[mbneu[i]]);
306	>	copycolor(bramp[i][1], mbRGB[mbneu[i]]);
307		}
308		/* compute color matrix */
309	<	for (i = 0; i < NMBMOD; i++)
310	<	for (j = 0; j < 3; j++) {
311	<	clrin[i][j] = bresp(picRGB[mbmod[i]][j]);
312	<	clrout[i][j] = mbRGB[mbmod[i]][j];
269	<	}
309	>	for (i = 0; i < NMBMOD; i++) {
310	>	bresp(clrin[i], picRGB[mbmod[i]]);
311	>	copycolor(clrout[i], mbRGB[mbmod[i]]);
312	>	}
313		compsoln(clrin, clrout, NMBMOD);
314	+	}
315	+
316	+
317	+	putmapping()                    /* put out color mapping for pcomb -f */
318	+	{
319	+	static char     cchar[3] = {'r', 'g', 'b'};
320	+	register int    i, j;
321		/* print brightness mapping */
322	<	printf("xval(i) : select(i, 0");
323	<	for (i = 0; i < NMBNEU; i++)
324	<	printf(", %g", bramp[i][0]);
325	<	printf(");\n");
326	<	printf("yval(i) : select(i, 0");
327	<	for (i = 0; i < NMBNEU; i++)
328	<	printf(", %g", bramp[i][1]);
329	<	printf(");\n");
330	<	printf("ifind(x,f,n) : if(1.5-n, 1, if(x-f(n), n, ifind(x,f,n-1)));\n");
331	<	printf("binterp(i,x) : ((xval(i+1)-x)*yval(i) + (x-xval(i))*yval(i+1))/\n");
332	<	printf("\t\t(xval(i+1) - xval(i));\n");
333	<	printf("bresp(x) : binterp(ifind(x,xval,xval(0)-1), x);\n");
334	<	printf("nred = bresp(ri(1));\n");
335	<	printf("ngrn = bresp(gi(1));\n");
336	<	printf("nblu = bresp(bi(1));\n");
322	>	for (j = 0; j < 3; j++) {
323	>	printf("%cxa(i) : select(i", cchar[j]);
324	>	for (i = 0; i < NMBNEU; i++)
325	>	printf(",%g", colval(bramp[i][0],j));
326	>	printf(");\n");
327	>	printf("%cya(i) : select(i", cchar[j]);
328	>	for (i = 0; i < NMBNEU; i++)
329	>	printf(",%g", colval(bramp[i][1],j));
330	>	printf(");\n");
331	>	printf("%c = %ci(1);\n", cchar[j], cchar[j]);
332	>	printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
333	>	cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
334	>	cchar[j], cchar[j]);
335	>	printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
336	>	printf("%cn = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
337	>	cchar[j], cchar[j], cchar[j],
338	>	cchar[j], cchar[j], cchar[j]);
339	>	printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
340	>	cchar[j], cchar[j], cchar[j],
341	>	cchar[j], cchar[j]);
342	>	printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
343	>	cchar[j], cchar[j], cchar[j], cchar[j]);
344	>	}
345		/* print color mapping */
346	<	printf("ro = %g*nred + %g*ngrn + %g*nblu\n",
347	<	solmat[0][0], solmat[1][0], solmat[2][0]);
348	<	printf("go = %g*nred + %g*ngrn + %g*nblu\n",
349	<	solmat[0][1], solmat[1][1], solmat[2][1]);
350	<	printf("bo = %g*nred + %g*ngrn + %g*nblu\n",
351	<	solmat[0][2], solmat[1][2], solmat[2][2]);
346	>	printf("ro = %g*rn + %g*gn + %g*bn ;\n",
347	>	solmat[0][0], solmat[0][1], solmat[0][2]);
348	>	printf("go = %g*rn + %g*gn + %g*bn ;\n",
349	>	solmat[1][0], solmat[1][1], solmat[1][2]);
350	>	printf("bo = %g*rn + %g*gn + %g*bn ;\n",
351	>	solmat[2][0], solmat[2][1], solmat[2][2]);
352		}
353
354
355	<	compsoln(cin, cout, n)          /* solve 3x3 system */
356	<	float   cin[][3], cout[][3];
355	>	compsoln(cin, cout, n)          /* solve 3xN system using least-squares */
356	>	COLOR   cin[], cout[];
357		int     n;
358		{
359		extern double   mx3d_adjoint(), fabs();
360		double  mat[3][3], invmat[3][3];
361		double  det;
362		double  colv[3], rowv[3];
363	<	register int    i, j;
363	>	register int    i, j, k;
364
365	<	if (n != 3) {
365	>	if (n < 3 | n > NMBMOD) {
366		fprintf(stderr, "%s: inconsistent code!\n", progname);
367		exit(1);
368		}
369	<	for (i = 0; i < 3; i++)
370	<	for (j = 0; j < 3; j++)
371	<	mat[i][j] = cin[j][i];
369	>	if (n == 3)
370	>	for (i = 0; i < 3; i++)
371	>	for (j = 0; j < 3; j++)
372	>	mat[i][j] = colval(cin[j],i);
373	>	else {                          /* compute A^t A */
374	>	for (i = 0; i < 3; i++)
375	>	for (j = i; j < 3; j++) {
376	>	mat[i][j] = 0.;
377	>	for (k = 0; k < n; k++)
378	>	mat[i][j] += colval(cin[k],i) *
379	>	colval(cin[k],j);
380	>	}
381	>	for (i = 1; i < 3; i++)         /* using symmetry */
382	>	for (j = 0; j < i; j++)
383	>	mat[i][j] = mat[j][i];
384	>	}
385		det = mx3d_adjoint(mat, invmat);
386		if (fabs(det) < 1e-4) {
387		fprintf(stderr, "%s: cannot compute color mapping\n",
#	Line 324 | Line 395 | int    n;
395		for (j = 0; j < 3; j++)
396		invmat[i][j] /= det;
397		for (i = 0; i < 3; i++) {
398	+	if (n == 3)
399	+	for (j = 0; j < 3; j++)
400	+	colv[j] = colval(cout[j],i);
401	+	else
402	+	for (j = 0; j < 3; j++) {
403	+	colv[j] = 0.;
404	+	for (k = 0; k < n; k++)
405	+	colv[j] += colval(cout[k],i) *
406	+	colval(cin[k],j);
407	+	}
408	+	mx3d_transform(colv, invmat, rowv);
409		for (j = 0; j < 3; j++)
410	<	rowv[j] = cout[j][i];
329	<	mx3d_transform(rowv, invmat, colv);
330	<	for (j = 0; j < 3; j++)
331	<	solmat[j][i] = colv[j];
410	>	solmat[i][j] = rowv[j];
411		}
412		}
413
#	Line 337 | Line 416 | cvtcolor(cout, cin)            /* convert color according to our
416		COLOR   cout, cin;
417		{
418		double  r, g, b;
340	–	double  r1, g1, b1;
419
420	<	r = bresp(colval(cin,RED));
421	<	g = bresp(colval(cin,GRN));
422	<	b = bresp(colval(cin,BLU));
423	<	r1 = r*solmat[0][0] + g*solmat[1][0] + b*solmat[2][0];
424	<	if (r1 < 0) r1 = 0;
425	<	g1 = r*solmat[0][1] + g*solmat[1][1] + b*solmat[2][1];
426	<	if (g1 < 0) g1 = 0;
427	<	b1 = r*solmat[0][2] + g*solmat[1][2] + b*solmat[2][2];
428	<	if (b1 < 0) b1 = 0;
429	<	setcolor(cout, r1, g1, b1);
420	>	bresp(cout, cin);
421	>	r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1]
422	>	+ colval(cout,2)*solmat[0][2];
423	>	if (r < 0) r = 0;
424	>	g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1]
425	>	+ colval(cout,2)*solmat[1][2];
426	>	if (g < 0) g = 0;
427	>	b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
428	>	+ colval(cout,2)*solmat[2][2];
429	>	if (b < 0) b = 0;
430	>	setcolor(cout, r, g, b);
431		}
432
433
#	Line 373 | Line 452 | putdebug()                     /* put out debugging picture */
452		exit(1);
453		}
454		/* finish debug header */
455	+	fputformat(COLRFMT, debugfp);
456		putc('\n', debugfp);
457		fprtresolu(xmax, ymax, debugfp);
458		for (y = ymax-1; y >= 0; y--) {

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)