src/px/macbethcal.c

/* Copyright (c) 1995 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Calibrate a scanned MacBeth Color Checker Chart
 *
 * Produce a .cal file suitable for use with pcomb.
 */

#include <stdio.h>
#ifdef MSDOS
#include <fcntl.h>
#endif
#include "color.h"
#include "resolu.h"
#include "pmap.h"

                                /* MacBeth colors */
#define DarkSkin        0
#define LightSkin       1
#define BlueSky         2
#define Foliage         3
#define BlueFlower      4
#define BluishGreen     5
#define Orange          6
#define PurplishBlue    7
#define ModerateRed     8
#define Purple          9
#define YellowGreen     10
#define OrangeYellow    11
#define Blue            12
#define Green           13
#define Red             14
#define Yellow          15
#define Magenta         16
#define Cyan            17
#define White           18
#define Neutral8        19
#define Neutral65       20
#define Neutral5        21
#define Neutral35       22
#define Black           23
                                /* computed from 5nm spectral measurements */
                                /* CIE 1931 2 degree obs, equal-energy white */
float   mbxyY[24][3] = {
                {0.462, 0.3769, 0.0932961},     /* DarkSkin */
                {0.4108, 0.3542, 0.410348},     /* LightSkin */
                {0.2626, 0.267, 0.181554},      /* BlueSky */
                {0.36, 0.4689, 0.108447},       /* Foliage */
                {0.2977, 0.2602, 0.248407},     /* BlueFlower */
                {0.2719, 0.3485, 0.401156},     /* BluishGreen */
                {0.52, 0.4197, 0.357899},       /* Orange */
                {0.229, 0.1866, 0.103911},      /* PurplishBlue */
                {0.4909, 0.3262, 0.242615},     /* ModerateRed */
                {0.3361, 0.2249, 0.0600102},    /* Purple */
                {0.3855, 0.4874, 0.42963},      /* YellowGreen */
                {0.4853, 0.4457, 0.476343},     /* OrangeYellow */
                {0.2026, 0.1369, 0.0529249},    /* Blue */
                {0.3007, 0.4822, 0.221226},     /* Green */
                {0.5805, 0.3238, 0.162167},     /* Red */
                {0.4617, 0.472, 0.64909},       /* Yellow */
                {0.4178, 0.2625, 0.233662},     /* Magenta */
                {0.2038, 0.2508, 0.167275},     /* Cyan */
                {0.3358, 0.337, 0.916877},      /* White */
                {0.3338, 0.3348, 0.604678},     /* Neutral.8 */
                {0.3333, 0.3349, 0.364566},     /* Neutral.65 */
                {0.3353, 0.3359, 0.200238},     /* Neutral.5 */
                {0.3363, 0.336, 0.0878721},     /* Neutral.35 */
                {0.3346, 0.3349, 0.0308383}     /* Black */
        };

COLOR   mbRGB[24];              /* MacBeth RGB values */

#define NMBNEU          6       /* Number of MacBeth neutral colors */
short   mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};

#define NMBMOD          17      /* Number of MacBeth unsaturated colors */
short   mbmod[NMBMOD] = {
                DarkSkin,LightSkin,BlueSky,Foliage,BlueFlower,BluishGreen,
                Orange,PurplishBlue,ModerateRed,Purple,YellowGreen,
                Black,Neutral35,Neutral5,Neutral65,Neutral8,White
        };

#define NMBSAT          8       /* Number of MacBeth saturated colors */
short   mbsat[NMBSAT] = {
                Red,Green,Blue,Magenta,Yellow,Cyan,
                Orange,OrangeYellow
        };

int     xmax, ymax;             /* input image dimensions */
int     bounds[4][2];           /* image coordinates of chart corners */
double  imgxfm[3][3];           /* coordinate transformation matrix */

COLOR   picRGB[24];             /* picture colors */

COLOR   bramp[NMBNEU][2];       /* brightness ramp (per primary) */
double  solmat[3][3];           /* color mapping matrix */

FILE    *debugfp = NULL;        /* debug output picture */
char    *progname;

extern char     *malloc();


main(argc, argv)
int     argc;
char    **argv;
{
        int     i;

        progname = argv[0];
        if (argc > 2 && !strcmp(argv[1], "-d")) {       /* debug output */
                if ((debugfp = fopen(argv[2], "w")) == NULL) {
                        perror(argv[2]);
                        exit(1);
                }
#ifdef MSDOS
                setmode(fileno(debugfp), O_BINARY);
#endif
                newheader("RADIANCE", debugfp);
                printargs(argc, argv, debugfp);
                argv += 2;
                argc -= 2;
        }
        if (argc != 3 && argc != 11)
                goto userr;
        if (strcmp(argv[1], "-") && freopen(argv[1], "r", stdin) == NULL) {
                perror(argv[1]);
                exit(1);
        }
        if (strcmp(argv[2], "-") && freopen(argv[2], "w", stdout) == NULL) {
                perror(argv[2]);
                exit(1);
        }
#ifdef MSDOS
        setmode(fileno(stdin), O_BINARY);
#endif
        if (checkheader(stdin, COLRFMT, NULL) < 0 ||
                        fgetresolu(&xmax, &ymax, stdin) < 0) {
                fprintf(stderr, "%s: bad input picture\n", progname);
                exit(1);
        }
                                        /* get chart boundaries */
        if (argc == 11) {
                for (i = 0; i < 4; i++) {
                        if (!isint(argv[2*i+3]) | !isint(argv[2*i+4]))
                                goto userr;
                        bounds[i][0] = atoi(argv[2*i+3]);
                        bounds[i][1] = atoi(argv[2*i+4]);
                }
        } else {
                bounds[0][0] = bounds[2][0] = .029*xmax + .5;
                bounds[0][1] = bounds[1][1] = .956*ymax + .5;
                bounds[1][0] = bounds[3][0] = .971*xmax + .5;
                bounds[2][1] = bounds[3][1] = .056*ymax + .5;
        }
        init();                         /* initialize */
        getcolors();                    /* get picture colors */
        compute();                      /* compute color mapping */
                                        /* print comment */
        printf("{ Color correction file computed by %s }\n", progname);
        printf("{ from scanned MacBetch color chart %s }\n", argv[1]);
        putmapping();                   /* put out color mapping */
        putdebug();                     /* put out debug picture */
        exit(0);
userr:
        fprintf(stderr, "Usage: %s [-d dbg.pic] input.pic output.cal [xul yul xur yur xll yll xlr ylr]\n",
                        progname);
        exit(1);
}


init()                          /* initialize */
{
        double  quad[4][2];
                                        /* make coordinate transformation */
        quad[0][0] = bounds[0][0];
        quad[0][1] = bounds[0][1];
        quad[1][0] = bounds[1][0];
        quad[1][1] = bounds[1][1];
        quad[2][0] = bounds[3][0];
        quad[2][1] = bounds[3][1];
        quad[3][0] = bounds[2][0];
        quad[3][1] = bounds[2][1];

        if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
                fprintf(stderr, "%s: bad chart boundaries\n", progname);
                exit(1);
        }
}


int
chartndx(x, y)                          /* find color number for position */
int     x, y;
{
        double  ipos[3], cpos[3];
        int     ix, iy;
        double  fx, fy;

        ipos[0] = x;
        ipos[1] = y;
        ipos[2] = 1;
        mx3d_transform(ipos, imgxfm, cpos);
        cpos[0] /= cpos[2];
        cpos[1] /= cpos[2];
        if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
                return(-1);
        ix = cpos[0];
        iy = cpos[1];
        fx = cpos[0] - ix;
        fy = cpos[1] - iy;
        if (fx < .35 || fx >= .65 || fy < .35 || fy >= .65)
                return(-1);
        return(iy*6 + ix);
}


getcolors()                             /* load in picture colors */
{
        COLR    *scanln;
        COLOR   pval;
        int     ccount[24];
        double  d;
        int     y;
        register int    x, i;

        scanln = (COLR *)malloc(xmax*sizeof(COLR));
        if (scanln == NULL) {
                perror(progname);
                exit(1);
        }
        for (i = 0; i < 24; i++) {
                setcolor(picRGB[i], 0., 0., 0.);
                ccount[i] = 0;
        }
        for (y = ymax-1; y >= 0; y--) {
                if (freadcolrs(scanln, xmax, stdin) < 0) {
                        fprintf(stderr, "%s: error reading input picture\n",
                                        progname);
                        exit(1);
                }
                for (x = 0; x < xmax; x++) {
                        i = chartndx(x, y);
                        if (i >= 0) {
                                colr_color(pval, scanln[x]);
                                addcolor(picRGB[i], pval);
                                ccount[i]++;
                        }
                }
        }
        for (i = 0; i < 24; i++) {
                if (ccount[i] == 0) {
                        fprintf(stderr, "%s: bad chart boundaries\n",
                                        progname);
                        exit(1);
                }
                d = 1.0/ccount[i];
                scalecolor(picRGB[i], d);
        }
        free((char *)scanln);
}


bresp(y, x)             /* piecewise linear interpolation of primaries */
COLOR   y, x;
{
        double  cv[3];
        register int    i, n;

        for (i = 0; i < 3; i++) {
                n = NMBNEU;
                while (n > 0 && colval(x,i) < colval(bramp[--n][0],i))
                        ;
                cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
                                                colval(bramp[n][1],i) +
                                (colval(x,i) - colval(bramp[n][0],i)) *
                                                colval(bramp[n+1][1],i)) /
                        (colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
                if (cv[i] < 0.) cv[i] = 0.;
        }
        setcolor(y, cv[0], cv[1], cv[2]);
}


compute()                       /* compute color mapping */
{
        COLOR   clrin[NMBMOD], clrout[NMBMOD];
        COLOR   ctmp;
        double  d;
        register int    i;
                                        /* map MacBeth colors to RGB space */
        for (i = 0; i < 24; i++) {
                d = mbxyY[i][2] / mbxyY[i][1];
                ctmp[0] = mbxyY[i][0] * d;
                ctmp[1] = mbxyY[i][2];
                ctmp[2] = (1. - mbxyY[i][0] - mbxyY[i][1]) * d;
                cie_rgb(mbRGB[i], ctmp);
        }
                                        /* compute piecewise luminance curve */
        for (i = 0; i < NMBNEU; i++) {
                copycolor(bramp[i][0], picRGB[mbneu[i]]);
                copycolor(bramp[i][1], mbRGB[mbneu[i]]);
        }
                                        /* compute color matrix */
        for (i = 0; i < NMBMOD; i++) {
                bresp(clrin[i], picRGB[mbmod[i]]);
                copycolor(clrout[i], mbRGB[mbmod[i]]);
        }
        compsoln(clrin, clrout, NMBMOD);
}


putmapping()                    /* put out color mapping for pcomb -f */
{
        static char     cchar[3] = {'r', 'g', 'b'};
        register int    i, j;
                                        /* print brightness mapping */
        for (j = 0; j < 3; j++) {
                printf("%cxa(i) : select(i", cchar[j]);
                for (i = 0; i < NMBNEU; i++)
                        printf(",%g", colval(bramp[i][0],j));
                printf(");\n");
                printf("%cya(i) : select(i", cchar[j]);
                for (i = 0; i < NMBNEU; i++)
                        printf(",%g", colval(bramp[i][1],j));
                printf(");\n");
                printf("%c = %ci(1);\n", cchar[j], cchar[j]);
                printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
                                cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
                                cchar[j], cchar[j]);
                printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
                printf("%cn = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
                                cchar[j], cchar[j], cchar[j],
                                cchar[j], cchar[j], cchar[j]);
                printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
                                cchar[j], cchar[j], cchar[j],
                                cchar[j], cchar[j]);
                printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
                                cchar[j], cchar[j], cchar[j], cchar[j]);
        }
                                        /* print color mapping */
        printf("ro = %g*rn + %g*gn + %g*bn ;\n",
                        solmat[0][0], solmat[0][1], solmat[0][2]);
        printf("go = %g*rn + %g*gn + %g*bn ;\n",
                        solmat[1][0], solmat[1][1], solmat[1][2]);
        printf("bo = %g*rn + %g*gn + %g*bn ;\n",
                        solmat[2][0], solmat[2][1], solmat[2][2]);
}


compsoln(cin, cout, n)          /* solve 3xN system using least-squares */
COLOR   cin[], cout[];
int     n;
{
        extern double   mx3d_adjoint(), fabs();
        double  mat[3][3], invmat[3][3];
        double  det;
        double  colv[3], rowv[3];
        register int    i, j, k;

        if (n < 3 | n > NMBMOD) {
                fprintf(stderr, "%s: inconsistent code!\n", progname);
                exit(1);
        }
        if (n == 3)
                for (i = 0; i < 3; i++)
                        for (j = 0; j < 3; j++)
                                mat[i][j] = colval(cin[j],i);
        else {                          /* compute A^t A */
                for (i = 0; i < 3; i++)
                        for (j = i; j < 3; j++) {
                                mat[i][j] = 0.;
                                for (k = 0; k < n; k++)
                                        mat[i][j] += colval(cin[k],i) *
                                                        colval(cin[k],j);
                        }
                for (i = 1; i < 3; i++)         /* using symmetry */
                        for (j = 0; j < i; j++)
                                mat[i][j] = mat[j][i];
        }
        det = mx3d_adjoint(mat, invmat);
        if (fabs(det) < 1e-4) {
                fprintf(stderr, "%s: cannot compute color mapping\n",
                                progname);
                solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
                solmat[0][1] = solmat[0][2] = solmat[1][0] =
                solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
                return;
        }
        for (i = 0; i < 3; i++)
                for (j = 0; j < 3; j++)
                        invmat[i][j] /= det;
        for (i = 0; i < 3; i++) {
                if (n == 3)
                        for (j = 0; j < 3; j++)
                                colv[j] = colval(cout[j],i);
                else
                        for (j = 0; j < 3; j++) {
                                colv[j] = 0.;
                                for (k = 0; k < n; k++)
                                        colv[j] += colval(cout[k],i) *
                                                        colval(cin[k],j);
                        }
                mx3d_transform(colv, invmat, rowv);
                for (j = 0; j < 3; j++)
                        solmat[i][j] = rowv[j];
        }
}


cvtcolor(cout, cin)             /* convert color according to our mapping */
COLOR   cout, cin;
{
        double  r, g, b;

        bresp(cout, cin);
        r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1]
                        + colval(cout,2)*solmat[0][2];
        if (r < 0) r = 0;
        g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1]
                        + colval(cout,2)*solmat[1][2];
        if (g < 0) g = 0;
        b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
                        + colval(cout,2)*solmat[2][2];
        if (b < 0) b = 0;
        setcolor(cout, r, g, b);
}


putdebug()                      /* put out debugging picture */
{
        COLOR   *scan;
        int     y;
        register int    x, i;

        if (debugfp == NULL)
                return;
        if (fseek(stdin, 0L, 0) == EOF) {
                fprintf(stderr, "%s: cannot seek on input picture\n", progname);
                exit(1);
        }
        getheader(stdin, NULL, NULL);           /* skip input header */
        fgetresolu(&xmax, &ymax, stdin);
                                                /* allocate scanline */
        scan = (COLOR *)malloc(xmax*sizeof(COLOR));
        if (scan == NULL) {
                perror(progname);
                exit(1);
        }
                                                /* finish debug header */
        putc('\n', debugfp);
        fprtresolu(xmax, ymax, debugfp);
        for (y = ymax-1; y >= 0; y--) {
                if (freadscan(scan, xmax, stdin) < 0) {
                        fprintf(stderr, "%s: error rereading input picture\n",
                                        progname);
                        exit(1);
                }
                for (x = 0; x < xmax; x++) {
                        i = chartndx(x, y);
                        if (i < 0)
                                cvtcolor(scan[x], scan[x]);
                        else
                                copycolor(scan[x], mbRGB[i]);
                }
                if (fwritescan(scan, xmax, debugfp) < 0) {
                        fprintf(stderr, "%s: error writing debugging picture\n",
                                        progname);
                        exit(1);
                }
        }
        free((char *)scan);
}
Revision:	2.4
Committed:	Wed Oct 11 20:48:17 1995 UTC (28 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.3:	+74 -62 lines
Log Message:	working version using least-squares matrix computation
#	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.4	#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	greg	2.1	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	greg	2.2	COLOR bramp[NMBNEU][2]; /* brightness ramp (per primary) */
100	greg	2.1	double solmat[3][3]; /* color mapping matrix */
101
102	greg	2.4	FILE debugfp = NULL; / debug output picture */
103	greg	2.1	char *progname;
104
105			extern char *malloc();
106
107
108			main(argc, argv)
109			int argc;
110			char **argv;
111			{
112			int i;
113
114			progname = argv[0];
115			if (argc > 2 && !strcmp(argv[1], "-d")) { /* debug output */
116			if ((debugfp = fopen(argv[2], "w")) == NULL) {
117			perror(argv[2]);
118			exit(1);
119			}
120			#ifdef MSDOS
121			setmode(fileno(debugfp), O_BINARY);
122			#endif
123			newheader("RADIANCE", debugfp);
124			printargs(argc, argv, debugfp);
125			argv += 2;
126			argc -= 2;
127			}
128			if (argc != 3 && argc != 11)
129			goto userr;
130			if (strcmp(argv[1], "-") && freopen(argv[1], "r", stdin) == NULL) {
131			perror(argv[1]);
132			exit(1);
133			}
134			if (strcmp(argv[2], "-") && freopen(argv[2], "w", stdout) == NULL) {
135			perror(argv[2]);
136			exit(1);
137			}
138			#ifdef MSDOS
139			setmode(fileno(stdin), O_BINARY);
140			#endif
141			if (checkheader(stdin, COLRFMT, NULL) < 0 \|\|
142			fgetresolu(&xmax, &ymax, stdin) < 0) {
143			fprintf(stderr, "%s: bad input picture\n", progname);
144			exit(1);
145			}
146			/* get chart boundaries */
147			if (argc == 11) {
148			for (i = 0; i < 4; i++) {
149			if (!isint(argv[2i+3]) \| !isint(argv[2i+4]))
150			goto userr;
151			bounds[i][0] = atoi(argv[2*i+3]);
152			bounds[i][1] = atoi(argv[2*i+4]);
153			}
154			} else {
155			bounds[0][0] = bounds[2][0] = .029*xmax + .5;
156			bounds[0][1] = bounds[1][1] = .956*ymax + .5;
157			bounds[1][0] = bounds[3][0] = .971*xmax + .5;
158			bounds[2][1] = bounds[3][1] = .056*ymax + .5;
159			}
160			init(); /* initialize */
161			getcolors(); /* get picture colors */
162	greg	2.2	compute(); /* compute color mapping */
163	greg	2.4	/* print comment */
164			printf("{ Color correction file computed by %s }\n", progname);
165			printf("{ from scanned MacBetch color chart %s }\n", argv[1]);
166	greg	2.1	putmapping(); /* put out color mapping */
167			putdebug(); /* put out debug picture */
168			exit(0);
169			userr:
170			fprintf(stderr, "Usage: %s [-d dbg.pic] input.pic output.cal [xul yul xur yur xll yll xlr ylr]\n",
171			progname);
172			exit(1);
173			}
174
175
176			init() /* initialize */
177			{
178			double quad[4][2];
179			/* make coordinate transformation */
180			quad[0][0] = bounds[0][0];
181			quad[0][1] = bounds[0][1];
182			quad[1][0] = bounds[1][0];
183			quad[1][1] = bounds[1][1];
184			quad[2][0] = bounds[3][0];
185			quad[2][1] = bounds[3][1];
186			quad[3][0] = bounds[2][0];
187			quad[3][1] = bounds[2][1];
188
189			if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
190			fprintf(stderr, "%s: bad chart boundaries\n", progname);
191			exit(1);
192			}
193			}
194
195
196			int
197			chartndx(x, y) /* find color number for position */
198			int x, y;
199			{
200			double ipos[3], cpos[3];
201			int ix, iy;
202			double fx, fy;
203
204			ipos[0] = x;
205			ipos[1] = y;
206			ipos[2] = 1;
207			mx3d_transform(ipos, imgxfm, cpos);
208			cpos[0] /= cpos[2];
209			cpos[1] /= cpos[2];
210			if (cpos[0] < 0. \|\| cpos[0] >= 6. \|\| cpos[1] < 0. \|\| cpos[1] >= 4.)
211			return(-1);
212			ix = cpos[0];
213			iy = cpos[1];
214			fx = cpos[0] - ix;
215			fy = cpos[1] - iy;
216			if (fx < .35 \|\| fx >= .65 \|\| fy < .35 \|\| fy >= .65)
217			return(-1);
218			return(iy*6 + ix);
219			}
220
221
222			getcolors() /* load in picture colors */
223			{
224			COLR *scanln;
225			COLOR pval;
226			int ccount[24];
227			double d;
228			int y;
229			register int x, i;
230
231			scanln = (COLR )malloc(xmaxsizeof(COLR));
232			if (scanln == NULL) {
233			perror(progname);
234			exit(1);
235			}
236			for (i = 0; i < 24; i++) {
237			setcolor(picRGB[i], 0., 0., 0.);
238			ccount[i] = 0;
239			}
240			for (y = ymax-1; y >= 0; y--) {
241			if (freadcolrs(scanln, xmax, stdin) < 0) {
242			fprintf(stderr, "%s: error reading input picture\n",
243			progname);
244			exit(1);
245			}
246			for (x = 0; x < xmax; x++) {
247			i = chartndx(x, y);
248			if (i >= 0) {
249			colr_color(pval, scanln[x]);
250			addcolor(picRGB[i], pval);
251			ccount[i]++;
252			}
253			}
254			}
255			for (i = 0; i < 24; i++) {
256			if (ccount[i] == 0) {
257			fprintf(stderr, "%s: bad chart boundaries\n",
258			progname);
259			exit(1);
260			}
261			d = 1.0/ccount[i];
262			scalecolor(picRGB[i], d);
263			}
264			free((char *)scanln);
265			}
266
267
268	greg	2.2	bresp(y, x) /* piecewise linear interpolation of primaries */
269			COLOR y, x;
270	greg	2.1	{
271	greg	2.3	double cv[3];
272	greg	2.2	register int i, n;
273	greg	2.1
274	greg	2.2	for (i = 0; i < 3; i++) {
275			n = NMBNEU;
276			while (n > 0 && colval(x,i) < colval(bramp[--n][0],i))
277			;
278	greg	2.3	cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
279	greg	2.2	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	greg	2.3	if (cv[i] < 0.) cv[i] = 0.;
284	greg	2.2	}
285	greg	2.3	setcolor(y, cv[0], cv[1], cv[2]);
286	greg	2.1	}
287
288
289	greg	2.2	compute() /* compute color mapping */
290	greg	2.1	{
291	greg	2.2	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	greg	2.1	/* compute piecewise luminance curve */
304			for (i = 0; i < NMBNEU; i++) {
305	greg	2.2	copycolor(bramp[i][0], picRGB[mbneu[i]]);
306			copycolor(bramp[i][1], mbRGB[mbneu[i]]);
307	greg	2.1	}
308			/* compute color matrix */
309	greg	2.2	for (i = 0; i < NMBMOD; i++) {
310			bresp(clrin[i], picRGB[mbmod[i]]);
311			copycolor(clrout[i], mbRGB[mbmod[i]]);
312			}
313	greg	2.1	compsoln(clrin, clrout, NMBMOD);
314	greg	2.2	}
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	greg	2.1	/* print brightness mapping */
322	greg	2.2	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	greg	2.1	/* print color mapping */
346	greg	2.2	printf("ro = %grn + %ggn + %g*bn ;\n",
347	greg	2.3	solmat[0][0], solmat[0][1], solmat[0][2]);
348	greg	2.2	printf("go = %grn + %ggn + %g*bn ;\n",
349	greg	2.3	solmat[1][0], solmat[1][1], solmat[1][2]);
350	greg	2.2	printf("bo = %grn + %ggn + %g*bn ;\n",
351	greg	2.3	solmat[2][0], solmat[2][1], solmat[2][2]);
352	greg	2.1	}
353
354
355	greg	2.4	compsoln(cin, cout, n) /* solve 3xN system using least-squares */
356	greg	2.2	COLOR cin[], cout[];
357	greg	2.1	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	greg	2.4	register int i, j, k;
364	greg	2.1
365	greg	2.4	if (n < 3 \| n > NMBMOD) {
366	greg	2.1	fprintf(stderr, "%s: inconsistent code!\n", progname);
367			exit(1);
368			}
369	greg	2.4	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	greg	2.1	det = mx3d_adjoint(mat, invmat);
386			if (fabs(det) < 1e-4) {
387			fprintf(stderr, "%s: cannot compute color mapping\n",
388			progname);
389			solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
390			solmat[0][1] = solmat[0][2] = solmat[1][0] =
391			solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
392			return;
393			}
394			for (i = 0; i < 3; i++)
395			for (j = 0; j < 3; j++)
396			invmat[i][j] /= det;
397			for (i = 0; i < 3; i++) {
398	greg	2.4	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	greg	2.3	mx3d_transform(colv, invmat, rowv);
409	greg	2.1	for (j = 0; j < 3; j++)
410	greg	2.3	solmat[i][j] = rowv[j];
411	greg	2.1	}
412			}
413
414	greg	2.3
415	greg	2.1	cvtcolor(cout, cin) /* convert color according to our mapping */
416			COLOR cout, cin;
417			{
418			double r, g, b;
419
420	greg	2.2	bresp(cout, cin);
421	greg	2.3	r = colval(cout,0)solmat[0][0] + colval(cout,1)solmat[0][1]
422			+ colval(cout,2)*solmat[0][2];
423	greg	2.2	if (r < 0) r = 0;
424	greg	2.3	g = colval(cout,0)solmat[1][0] + colval(cout,1)solmat[1][1]
425			+ colval(cout,2)*solmat[1][2];
426	greg	2.2	if (g < 0) g = 0;
427	greg	2.3	b = colval(cout,0)solmat[2][0] + colval(cout,1)solmat[2][1]
428	greg	2.2	+ colval(cout,2)*solmat[2][2];
429			if (b < 0) b = 0;
430			setcolor(cout, r, g, b);
431	greg	2.1	}
432
433
434			putdebug() /* put out debugging picture */
435			{
436			COLOR *scan;
437			int y;
438			register int x, i;
439
440			if (debugfp == NULL)
441			return;
442			if (fseek(stdin, 0L, 0) == EOF) {
443			fprintf(stderr, "%s: cannot seek on input picture\n", progname);
444			exit(1);
445			}
446			getheader(stdin, NULL, NULL); /* skip input header */
447			fgetresolu(&xmax, &ymax, stdin);
448			/* allocate scanline */
449			scan = (COLOR )malloc(xmaxsizeof(COLOR));
450			if (scan == NULL) {
451			perror(progname);
452			exit(1);
453			}
454			/* finish debug header */
455			putc('\n', debugfp);
456			fprtresolu(xmax, ymax, debugfp);
457			for (y = ymax-1; y >= 0; y--) {
458			if (freadscan(scan, xmax, stdin) < 0) {
459			fprintf(stderr, "%s: error rereading input picture\n",
460			progname);
461			exit(1);
462			}
463			for (x = 0; x < xmax; x++) {
464			i = chartndx(x, y);
465			if (i < 0)
466			cvtcolor(scan[x], scan[x]);
467			else
468			copycolor(scan[x], mbRGB[i]);
469			}
470			if (fwritescan(scan, xmax, debugfp) < 0) {
471			fprintf(stderr, "%s: error writing debugging picture\n",
472			progname);
473			exit(1);
474			}
475			}
476			free((char *)scan);
477			}