Back to index

radiance  4R0+20100331
macbethcal.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: macbethcal.c,v 2.25 2008/11/10 19:08:19 greg Exp $";
00003 #endif
00004 /*
00005  * Calibrate a scanned MacBeth Color Checker Chart
00006  *
00007  * Produce a .cal file suitable for use with pcomb,
00008  * or .cwp file suitable for use with pcwarp.
00009  *
00010  * Warping code depends on conformance of COLOR and W3VEC types.
00011  */
00012 
00013 #include <stdio.h>
00014 #include <math.h>
00015 #include <time.h>
00016 
00017 #include "platform.h"
00018 #include "rtprocess.h"
00019 #include "rtio.h"
00020 #include "color.h"
00021 #include "resolu.h"
00022 #include "pmap.h"
00023 #include "warp3d.h"
00024 #include "mx3.h"
00025 
00026                             /* MacBeth colors */
00027 #define       DarkSkin      0
00028 #define       LightSkin     1
00029 #define       BlueSky              2
00030 #define       Foliage              3
00031 #define       BlueFlower    4
00032 #define       BluishGreen   5
00033 #define       Orange        6
00034 #define       PurplishBlue  7
00035 #define       ModerateRed   8
00036 #define       Purple        9
00037 #define       YellowGreen   10
00038 #define       OrangeYellow  11
00039 #define       Blue          12
00040 #define       Green         13
00041 #define       Red           14
00042 #define       Yellow        15
00043 #define       Magenta              16
00044 #define       Cyan          17
00045 #define       White         18
00046 #define       Neutral8      19
00047 #define       Neutral65     20
00048 #define       Neutral5      21
00049 #define       Neutral35     22
00050 #define       Black         23
00051                             /* computed from 5nm spectral measurements */
00052                             /* CIE 1931 2 degree obs, equal-energy white */
00053 float  mbxyY[24][3] = {
00054               {0.462, 0.3769, 0.0932961}, /* DarkSkin */
00055               {0.4108, 0.3542, 0.410348}, /* LightSkin */
00056               {0.2626, 0.267, 0.181554},  /* BlueSky */
00057               {0.36, 0.4689, 0.108447},   /* Foliage */
00058               {0.2977, 0.2602, 0.248407}, /* BlueFlower */
00059               {0.2719, 0.3485, 0.401156}, /* BluishGreen */
00060               {0.52, 0.4197, 0.357899},   /* Orange */
00061               {0.229, 0.1866, 0.103911},  /* PurplishBlue */
00062               {0.4909, 0.3262, 0.242615}, /* ModerateRed */
00063               {0.3361, 0.2249, 0.0600102},       /* Purple */
00064               {0.3855, 0.4874, 0.42963},  /* YellowGreen */
00065               {0.4853, 0.4457, 0.476343}, /* OrangeYellow */
00066               {0.2026, 0.1369, 0.0529249},       /* Blue */
00067               {0.3007, 0.4822, 0.221226}, /* Green */
00068               {0.5805, 0.3238, 0.162167}, /* Red */
00069               {0.4617, 0.472, 0.64909},   /* Yellow */
00070               {0.4178, 0.2625, 0.233662}, /* Magenta */
00071               {0.2038, 0.2508, 0.167275}, /* Cyan */
00072               {0.3358, 0.337, 0.916877},  /* White */
00073               {0.3338, 0.3348, 0.604678}, /* Neutral.8 */
00074               {0.3333, 0.3349, 0.364566}, /* Neutral.65 */
00075               {0.3353, 0.3359, 0.200238}, /* Neutral.5 */
00076               {0.3363, 0.336, 0.0878721}, /* Neutral.35 */
00077               {0.3346, 0.3349, 0.0308383} /* Black */
00078        };
00079 
00080 COLOR  mbRGB[24];           /* MacBeth RGB values */
00081 
00082 #define       NMBNEU        6      /* Number of MacBeth neutral colors */
00083 short  mbneu[NMBNEU] = {Black,Neutral35,Neutral5,Neutral65,Neutral8,White};
00084 
00085 #define  NEUFLGS     (1L<<White|1L<<Neutral8|1L<<Neutral65| \
00086                             1L<<Neutral5|1L<<Neutral35|1L<<Black)
00087 
00088 #define  SATFLGS     (1L<<Red|1L<<Green|1L<<Blue|1L<<Magenta|1L<<Yellow| \
00089                      1L<<Cyan|1L<<Orange|1L<<Purple|1L<<PurplishBlue| \
00090                      1L<<YellowGreen|1<<OrangeYellow|1L<<BlueFlower)
00091 
00092 #define  UNSFLGS     (1L<<DarkSkin|1L<<LightSkin|1L<<BlueSky|1L<<Foliage| \
00093                      1L<<BluishGreen|1L<<ModerateRed)
00094 
00095 #define  REQFLGS     NEUFLGS                     /* need these colors */
00096 #define  MODFLGS     (NEUFLGS|UNSFLGS)    /* should be in gamut */
00097 
00098 #define  RG_BORD     0      /* patch border */
00099 #define  RG_CENT     01     /* central region of patch */
00100 #define  RG_ORIG     02     /* original color region */
00101 #define  RG_CORR     04     /* corrected color region */
00102 
00103 #ifndef  DISPCOM
00104 #define  DISPCOM     "ximage -op \"%s\""
00105 #endif
00106 
00107 int    scanning = 1;        /* scanned input (or recorded output)? */
00108 double irrad = 1.0;         /* irradiance multiplication factor */
00109 int    rawmap = 0;          /* put out raw color mapping? */
00110 
00111 int    xmax, ymax;          /* input image dimensions */
00112 int    bounds[4][2];        /* image coordinates of chart corners */
00113 double imgxfm[3][3];        /* coordinate transformation matrix */
00114 
00115 COLOR  inpRGB[24];          /* measured or scanned input colors */
00116 long   inpflags = 0;        /* flags of which colors were input */
00117 long   gmtflags = 0;        /* flags of out-of-gamut colors */
00118 
00119 COLOR  bramp[NMBNEU][2];    /* brightness ramp (per primary) */
00120 COLORMAT      solmat;              /* color mapping matrix */
00121 COLOR  colmin, colmax;             /* gamut limits */
00122 
00123 WARP3D *wcor = NULL;        /* color space warp */
00124 
00125 FILE   *debugfp = NULL;     /* debug output picture */
00126 char   *progname;
00127 
00128 static void init(void);
00129 static int chartndx(int     x, int y, int *np);
00130 static void getpicture(void);
00131 static void getcolors(void);
00132 static void bresp(COLOR     y, COLOR      x);
00133 static void compute(void);
00134 static void putmapping(void);
00135 static void compsoln(COLOR  cin[], COLOR  cout[], int   n);
00136 static void cwarp(void);
00137 static int cvtcolor(COLOR   cout, COLOR   cin);
00138 static int cresp(COLOR      cout, COLOR   cin);
00139 static void xyY2RGB(COLOR   rgbout, float xyYin[3]);
00140 static void picdebug(void);
00141 static void clrdebug(void);
00142 static void getpos(char     *name, int    bnds[2], FILE *fp);
00143 static void pickchartpos(char      *pfn);
00144 
00145 
00146 int
00147 main(
00148        int    argc,
00149        char   **argv
00150 )
00151 {
00152        int    i;
00153 
00154        progname = argv[0];
00155        for (i = 1; i < argc && argv[i][0] == '-'; i++)
00156               switch (argv[i][1]) {
00157               case 'd':                          /* debug output */
00158                      i++;
00159                      if (badarg(argc-i, argv+i, "s"))
00160                             goto userr;
00161                      if ((debugfp = fopen(argv[i], "w")) == NULL) {
00162                             perror(argv[i]);
00163                             exit(1);
00164                      }
00165                      SET_FILE_BINARY(debugfp);
00166                      newheader("RADIANCE", debugfp);           /* start */
00167                      printargs(argc, argv, debugfp);           /* header */
00168                      break;
00169               case 'p':                          /* picture position */
00170                      if (badarg(argc-i-1, argv+i+1, "iiiiiiii"))
00171                             goto userr;
00172                      bounds[0][0] = atoi(argv[++i]);
00173                      bounds[0][1] = atoi(argv[++i]);
00174                      bounds[1][0] = atoi(argv[++i]);
00175                      bounds[1][1] = atoi(argv[++i]);
00176                      bounds[2][0] = atoi(argv[++i]);
00177                      bounds[2][1] = atoi(argv[++i]);
00178                      bounds[3][0] = atoi(argv[++i]);
00179                      bounds[3][1] = atoi(argv[++i]);
00180                      scanning = 2;
00181                      break;
00182               case 'P':                          /* pick position */
00183                      scanning = 3;
00184                      break;
00185               case 'i':                          /* irradiance factor */
00186                      i++;
00187                      if (badarg(argc-i, argv+i, "f"))
00188                             goto userr;
00189                      irrad = atof(argv[i]);
00190                      break;
00191               case 'm':                          /* raw map output */
00192                      rawmap = 1;
00193                      break;
00194               case 'c':                          /* color input */
00195                      scanning = 0;
00196                      break;
00197               default:
00198                      goto userr;
00199               }
00200                                                  /* open files */
00201        if (i < argc && freopen(argv[i], "r", stdin) == NULL) {
00202               perror(argv[i]);
00203               exit(1);
00204        }
00205        if (i+1 < argc && freopen(argv[i+1], "w", stdout) == NULL) {
00206               perror(argv[i+1]);
00207               exit(1);
00208        }
00209        if (scanning) {                    /* load input picture header */
00210               SET_FILE_BINARY(stdin);
00211               if (checkheader(stdin, COLRFMT, NULL) < 0 ||
00212                             fgetresolu(&xmax, &ymax, stdin) < 0) {
00213                      fprintf(stderr, "%s: bad input picture\n", progname);
00214                      exit(1);
00215               }
00216               if (scanning == 3) {
00217                      if (i >= argc)
00218                             goto userr;
00219                      pickchartpos(argv[i]);
00220                      scanning = 2;
00221               }
00222        } else {                    /* else set default xmax and ymax */
00223               xmax = 512;
00224               ymax = 2*512/3;
00225        }
00226        if (scanning != 2) {        /* use default boundaries */
00227               bounds[0][0] = bounds[2][0] = .029*xmax + .5;
00228               bounds[0][1] = bounds[1][1] = .956*ymax + .5;
00229               bounds[1][0] = bounds[3][0] = .971*xmax + .5;
00230               bounds[2][1] = bounds[3][1] = .056*ymax + .5;
00231        }
00232        init();                            /* initialize */
00233        if (scanning)               /* get picture colors */
00234               getpicture();
00235        else
00236               getcolors();
00237        compute();                  /* compute color mapping */
00238        if (rawmap) {               /* print out raw correspondence */
00239               register int  j;
00240 
00241               printf("# Color correspondence produced by:\n#\t\t");
00242               printargs(argc, argv, stdout);
00243               printf("#\tUsage: pcwarp %s uncorrected.hdr > corrected.hdr\n",
00244                             i+1 < argc ? argv[i+1] : "{this_file}");
00245               printf("#\t   Or: pcond [options] -m %s orig.hdr > output.hdr\n",
00246                             i+1 < argc ? argv[i+1] : "{this_file}");
00247               for (j = 0; j < 24; j++)
00248                      printf("%f %f %f    %f %f %f\n",
00249                             colval(inpRGB[j],RED), colval(inpRGB[j],GRN),
00250                             colval(inpRGB[j],BLU), colval(mbRGB[j],RED),
00251                             colval(mbRGB[j],GRN), colval(mbRGB[j],BLU));
00252               if (scanning && debugfp != NULL)
00253                      cwarp();             /* color warp for debugging */
00254        } else {                    /* print color mapping */
00255                                           /* print header */
00256               printf("{\n\tColor correction file computed by:\n\t\t");
00257               printargs(argc, argv, stdout);
00258               printf("\n\tUsage: pcomb -f %s uncorrected.hdr > corrected.hdr\n",
00259                             i+1 < argc ? argv[i+1] : "{this_file}");
00260               if (!scanning)
00261                      printf("\t   Or: pcond [options] -f %s orig.hdr > output.hdr\n",
00262                                    i+1 < argc ? argv[i+1] : "{this_file}");
00263               printf("}\n");
00264               putmapping();               /* put out color mapping */
00265        }
00266        if (debugfp != NULL) {             /* put out debug picture */
00267               if (scanning)
00268                      picdebug();
00269               else
00270                      clrdebug();
00271        }
00272        exit(0);
00273 userr:
00274        fprintf(stderr,
00275 "Usage: %s [-d dbg.hdr][-P | -p xul yul xur yur xll yll xlr ylr][-i irrad][-m] input.hdr [output.{cal|cwp}]\n",
00276                      progname);
00277        fprintf(stderr, "   or: %s [-d dbg.hdr][-i irrad][-m] -c [xyY.dat [output.{cal|cwp}]]\n",
00278                      progname);
00279        exit(1);
00280        return 1; /* pro forma return */
00281 }
00282 
00283 
00284 static void
00285 init(void)                         /* initialize */
00286 {
00287        double quad[4][2];
00288        register int  i;
00289                                    /* make coordinate transformation */
00290        quad[0][0] = bounds[0][0];
00291        quad[0][1] = bounds[0][1];
00292        quad[1][0] = bounds[1][0];
00293        quad[1][1] = bounds[1][1];
00294        quad[2][0] = bounds[3][0];
00295        quad[2][1] = bounds[3][1];
00296        quad[3][0] = bounds[2][0];
00297        quad[3][1] = bounds[2][1];
00298 
00299        if (pmap_quad_rect(0., 0., 6., 4., quad, imgxfm) == PMAP_BAD) {
00300               fprintf(stderr, "%s: bad chart boundaries\n", progname);
00301               exit(1);
00302        }
00303                                    /* map MacBeth colors to RGB space */
00304        for (i = 0; i < 24; i++) {
00305               xyY2RGB(mbRGB[i], mbxyY[i]);
00306               scalecolor(mbRGB[i], irrad);
00307        }
00308 }
00309 
00310 
00311 static int
00312 chartndx(                   /* find color number for position */
00313        int    x,
00314        int y,
00315        int    *np
00316 )
00317 {
00318        double ipos[3], cpos[3];
00319        int    ix, iy;
00320        double fx, fy;
00321 
00322        ipos[0] = x;
00323        ipos[1] = y;
00324        ipos[2] = 1;
00325        mx3d_transform(ipos, imgxfm, cpos);
00326        cpos[0] /= cpos[2];
00327        cpos[1] /= cpos[2];
00328        if (cpos[0] < 0. || cpos[0] >= 6. || cpos[1] < 0. || cpos[1] >= 4.)
00329               return(RG_BORD);
00330        ix = cpos[0];
00331        iy = cpos[1];
00332        fx = cpos[0] - ix;
00333        fy = cpos[1] - iy;
00334        *np = iy*6 + ix;
00335        if (fx >= 0.35 && fx < 0.65 && fy >= 0.35 && fy < 0.65)
00336               return(RG_CENT);
00337        if (fx < 0.05 || fx >= 0.95 || fy < 0.05 || fy >= 0.95)
00338               return(RG_BORD);
00339        if (fx >= 0.5)                     /* right side is corrected */
00340               return(RG_CORR);
00341        return(RG_ORIG);            /* left side is original */
00342 }
00343 
00344 
00345 static void
00346 getpicture(void)                          /* load in picture colors */
00347 {
00348        COLR   *scanln;
00349        COLOR  pval;
00350        int    ccount[24];
00351        double d;
00352        int    y, i;
00353        register int  x;
00354 
00355        scanln = (COLR *)malloc(xmax*sizeof(COLR));
00356        if (scanln == NULL) {
00357               perror(progname);
00358               exit(1);
00359        }
00360        for (i = 0; i < 24; i++) {
00361               setcolor(inpRGB[i], 0., 0., 0.);
00362               ccount[i] = 0;
00363        }
00364        for (y = ymax-1; y >= 0; y--) {
00365               if (freadcolrs(scanln, xmax, stdin) < 0) {
00366                      fprintf(stderr, "%s: error reading input picture\n",
00367                                    progname);
00368                      exit(1);
00369               }
00370               for (x = 0; x < xmax; x++)
00371                      if (chartndx(x, y, &i) == RG_CENT) {
00372                             colr_color(pval, scanln[x]);
00373                             addcolor(inpRGB[i], pval);
00374                             ccount[i]++;
00375                      }
00376        }
00377        for (i = 0; i < 24; i++) {         /* compute averages */
00378               if (ccount[i] == 0)
00379                      continue;
00380               d = 1./ccount[i];
00381               scalecolor(inpRGB[i], d);
00382               inpflags |= 1L<<i;
00383        }
00384        free((void *)scanln);
00385 }
00386 
00387 
00388 static void
00389 getcolors(void)                    /* get xyY colors from standard input */
00390 {
00391        int    gotwhite = 0;
00392        COLOR  whiteclr;
00393        int    n;
00394        float  xyYin[3];
00395 
00396        while (fgetval(stdin, 'i', &n) == 1) {           /* read colors */
00397               if ((n < 0) | (n > 24) ||
00398                             fgetval(stdin, 'f', &xyYin[0]) != 1 ||
00399                             fgetval(stdin, 'f', &xyYin[1]) != 1 ||
00400                             fgetval(stdin, 'f', &xyYin[2]) != 1 ||
00401                             (xyYin[0] < 0.) | (xyYin[1] < 0.) ||
00402                             xyYin[0] + xyYin[1] > 1.) {
00403                      fprintf(stderr, "%s: bad color input data\n",
00404                                    progname);
00405                      exit(1);
00406               }
00407               if (n == 0) {                      /* calibration white */
00408                      xyY2RGB(whiteclr, xyYin);
00409                      gotwhite++;
00410               } else {                           /* standard color */
00411                      n--;
00412                      xyY2RGB(inpRGB[n], xyYin);
00413                      inpflags |= 1L<<n;
00414               }
00415        }
00416                                    /* normalize colors */
00417        if (!gotwhite) {
00418               if (!(inpflags & 1L<<White)) {
00419                      fprintf(stderr, "%s: missing input for White\n",
00420                                    progname);
00421                      exit(1);
00422               }
00423               setcolor(whiteclr,
00424                      colval(inpRGB[White],RED)/colval(mbRGB[White],RED),
00425                      colval(inpRGB[White],GRN)/colval(mbRGB[White],GRN),
00426                      colval(inpRGB[White],BLU)/colval(mbRGB[White],BLU));
00427        }
00428        for (n = 0; n < 24; n++)
00429               if (inpflags & 1L<<n)
00430                      setcolor(inpRGB[n],
00431                             colval(inpRGB[n],RED)/colval(whiteclr,RED),
00432                             colval(inpRGB[n],GRN)/colval(whiteclr,GRN),
00433                             colval(inpRGB[n],BLU)/colval(whiteclr,BLU));
00434 }
00435 
00436 
00437 static void
00438 bresp(        /* piecewise linear interpolation of primaries */
00439        COLOR  y,
00440        COLOR  x
00441 )
00442 {
00443        register int  i, n;
00444 
00445        for (i = 0; i < 3; i++) {
00446               for (n = 0; n < NMBNEU-2; n++)
00447                      if (colval(x,i) < colval(bramp[n+1][0],i))
00448                             break;
00449               colval(y,i) = ((colval(bramp[n+1][0],i) - colval(x,i)) *
00450                                           colval(bramp[n][1],i) +
00451                             (colval(x,i) - colval(bramp[n][0],i)) *
00452                                           colval(bramp[n+1][1],i)) /
00453                      (colval(bramp[n+1][0],i) - colval(bramp[n][0],i));
00454        }
00455 }
00456 
00457 
00458 static void
00459 compute(void)               /* compute color mapping */
00460 {
00461        COLOR  clrin[24], clrout[24];
00462        long   cflags;
00463        COLOR  ctmp;
00464        register int  i, n;
00465                                    /* did we get what we need? */
00466        if ((inpflags & REQFLGS) != REQFLGS) {
00467               fprintf(stderr, "%s: missing required input colors\n",
00468                             progname);
00469               exit(1);
00470        }
00471                                    /* compute piecewise luminance curve */
00472        for (i = 0; i < NMBNEU; i++) {
00473               copycolor(bramp[i][0], inpRGB[mbneu[i]]);
00474               for (n = i ? 3 : 0; n--; )
00475                      if (colval(bramp[i][0],n) <=
00476                                    colval(bramp[i-1][0],n)+1e-7) {
00477                             fprintf(stderr,
00478               "%s: non-increasing neutral patch\n", progname);
00479                             exit(1);
00480                      }
00481               copycolor(bramp[i][1], mbRGB[mbneu[i]]);
00482        }
00483                                    /* compute color space gamut */
00484        if (scanning) {
00485               copycolor(colmin, cblack);
00486               copycolor(colmax, cwhite);
00487               scalecolor(colmax, irrad);
00488        } else
00489               for (i = 0; i < 3; i++) {
00490                      colval(colmin,i) = colval(bramp[0][0],i) -
00491                             colval(bramp[0][1],i) *
00492                             (colval(bramp[1][0],i)-colval(bramp[0][0],i)) /
00493                             (colval(bramp[1][1],i)-colval(bramp[1][0],i));
00494                      colval(colmax,i) = colval(bramp[NMBNEU-2][0],i) +
00495                             (1.-colval(bramp[NMBNEU-2][1],i)) *
00496                             (colval(bramp[NMBNEU-1][0],i) -
00497                                    colval(bramp[NMBNEU-2][0],i)) /
00498                             (colval(bramp[NMBNEU-1][1],i) -
00499                                    colval(bramp[NMBNEU-2][1],i));
00500               }
00501                                    /* compute color mapping */
00502        do {
00503               cflags = inpflags & ~gmtflags;
00504               n = 0;                      /* compute transform matrix */
00505               for (i = 0; i < 24; i++)
00506                      if (cflags & 1L<<i) {
00507                             bresp(clrin[n], inpRGB[i]);
00508                             copycolor(clrout[n], mbRGB[i]);
00509                             n++;
00510                      }
00511               compsoln(clrin, clrout, n);
00512               if (irrad > 0.99 && irrad < 1.01)  /* check gamut */
00513                      for (i = 0; i < 24; i++)
00514                             if (cflags & 1L<<i && cvtcolor(ctmp, mbRGB[i]))
00515                                    gmtflags |= 1L<<i;
00516        } while (cflags & gmtflags);
00517        if (gmtflags & MODFLGS)
00518               fprintf(stderr,
00519               "%s: warning - some moderate colors are out of gamut\n",
00520                             progname);
00521 }
00522 
00523 
00524 static void
00525 putmapping(void)                   /* put out color mapping */
00526 {
00527        static char   cchar[3] = {'r', 'g', 'b'};
00528        register int  i, j;
00529                                    /* print brightness mapping */
00530        for (j = 0; j < 3; j++) {
00531               printf("%cxa(i) : select(i", cchar[j]);
00532               for (i = 0; i < NMBNEU; i++)
00533                      printf(",%g", colval(bramp[i][0],j));
00534               printf(");\n");
00535               printf("%cya(i) : select(i", cchar[j]);
00536               for (i = 0; i < NMBNEU; i++)
00537                      printf(",%g", colval(bramp[i][1],j));
00538               printf(");\n");
00539               printf("%cfi(n) = if(n-%g, %d, if(%cxa(n+1)-%c, n, %cfi(n+1)));\n",
00540                             cchar[j], NMBNEU-1.5, NMBNEU-1, cchar[j],
00541                             cchar[j], cchar[j]);
00542               printf("%cndx = %cfi(1);\n", cchar[j], cchar[j]);
00543               printf("%c%c = ((%cxa(%cndx+1)-%c)*%cya(%cndx) + ",
00544                             cchar[j], scanning?'n':'o', cchar[j],
00545                             cchar[j], cchar[j], cchar[j], cchar[j]);
00546               printf("(%c-%cxa(%cndx))*%cya(%cndx+1)) /\n",
00547                             cchar[j], cchar[j], cchar[j],
00548                             cchar[j], cchar[j]);
00549               printf("\t\t(%cxa(%cndx+1) - %cxa(%cndx)) ;\n",
00550                             cchar[j], cchar[j], cchar[j], cchar[j]);
00551        }
00552                                    /* print color mapping */
00553        if (scanning) {
00554               printf("r = ri(1); g = gi(1); b = bi(1);\n");
00555               printf("ro = %g*rn + %g*gn + %g*bn ;\n",
00556                             solmat[0][0], solmat[0][1], solmat[0][2]);
00557               printf("go = %g*rn + %g*gn + %g*bn ;\n",
00558                             solmat[1][0], solmat[1][1], solmat[1][2]);
00559               printf("bo = %g*rn + %g*gn + %g*bn ;\n",
00560                             solmat[2][0], solmat[2][1], solmat[2][2]);
00561        } else {
00562               printf("r1 = ri(1); g1 = gi(1); b1 = bi(1);\n");
00563               printf("r = %g*r1 + %g*g1 + %g*b1 ;\n",
00564                             solmat[0][0], solmat[0][1], solmat[0][2]);
00565               printf("g = %g*r1 + %g*g1 + %g*b1 ;\n",
00566                             solmat[1][0], solmat[1][1], solmat[1][2]);
00567               printf("b = %g*r1 + %g*g1 + %g*b1 ;\n",
00568                             solmat[2][0], solmat[2][1], solmat[2][2]);
00569        }
00570 }
00571 
00572 
00573 static void
00574 compsoln(            /* solve 3xN system using least-squares */
00575        COLOR  cin[],
00576        COLOR  cout[],
00577        int    n
00578 )
00579 {
00580        extern double mx3d_adjoint(), fabs();
00581        double mat[3][3], invmat[3][3];
00582        double det;
00583        double colv[3], rowv[3];
00584        register int  i, j, k;
00585 
00586        if (n < 3) {
00587               fprintf(stderr, "%s: too few colors to match!\n", progname);
00588               exit(1);
00589        }
00590        if (n == 3)
00591               for (i = 0; i < 3; i++)
00592                      for (j = 0; j < 3; j++)
00593                             mat[i][j] = colval(cin[j],i);
00594        else {                      /* compute A^t A */
00595               for (i = 0; i < 3; i++)
00596                      for (j = i; j < 3; j++) {
00597                             mat[i][j] = 0.;
00598                             for (k = 0; k < n; k++)
00599                                    mat[i][j] += colval(cin[k],i) *
00600                                                  colval(cin[k],j);
00601                      }
00602               for (i = 1; i < 3; i++)            /* using symmetry */
00603                      for (j = 0; j < i; j++)
00604                             mat[i][j] = mat[j][i];
00605        }
00606        det = mx3d_adjoint(mat, invmat);
00607        if (fabs(det) < 1e-4) {
00608               fprintf(stderr, "%s: cannot compute color mapping\n",
00609                             progname);
00610               solmat[0][0] = solmat[1][1] = solmat[2][2] = 1.;
00611               solmat[0][1] = solmat[0][2] = solmat[1][0] =
00612               solmat[1][2] = solmat[2][0] = solmat[2][1] = 0.;
00613               return;
00614        }
00615        for (i = 0; i < 3; i++)
00616               for (j = 0; j < 3; j++)
00617                      invmat[i][j] /= det;
00618        for (i = 0; i < 3; i++) {
00619               if (n == 3)
00620                      for (j = 0; j < 3; j++)
00621                             colv[j] = colval(cout[j],i);
00622               else
00623                      for (j = 0; j < 3; j++) {
00624                             colv[j] = 0.;
00625                             for (k = 0; k < n; k++)
00626                                    colv[j] += colval(cout[k],i) *
00627                                                  colval(cin[k],j);
00628                      }
00629               mx3d_transform(colv, invmat, rowv);
00630               for (j = 0; j < 3; j++)
00631                      solmat[i][j] = rowv[j];
00632        }
00633 }
00634 
00635 
00636 static void
00637 cwarp(void)                        /* compute color warp map */
00638 {
00639        register int  i;
00640 
00641        if ((wcor = new3dw(W3EXACT)) == NULL)
00642               goto memerr;
00643        for (i = 0; i < 24; i++)
00644               if (!add3dpt(wcor, inpRGB[i], mbRGB[i]))
00645                      goto memerr;
00646        return;
00647 memerr:
00648        perror(progname);
00649        exit(1);
00650 }
00651 
00652 
00653 static int
00654 cvtcolor(            /* convert color according to our mapping */
00655        COLOR  cout,
00656        COLOR  cin
00657 )
00658 {
00659        COLOR  ctmp;
00660        int    clipped;
00661 
00662        if (wcor != NULL) {
00663               clipped = warp3d(cout, cin, wcor);
00664               clipped |= clipgamut(cout,bright(cout),CGAMUT,colmin,colmax);
00665        } else if (scanning) {
00666               bresp(ctmp, cin);
00667               clipped = cresp(cout, ctmp);
00668        } else {
00669               clipped = cresp(ctmp, cin);
00670               bresp(cout, ctmp);
00671        }
00672        return(clipped);
00673 }
00674 
00675 
00676 static int
00677 cresp(        /* transform color according to matrix */
00678        COLOR  cout,
00679        COLOR  cin
00680 )
00681 {
00682        colortrans(cout, solmat, cin);
00683        return(clipgamut(cout, bright(cout), CGAMUT, colmin, colmax));
00684 }
00685 
00686 
00687 static void
00688 xyY2RGB(             /* convert xyY to RGB */
00689        COLOR  rgbout,
00690        register float       xyYin[3]
00691 )
00692 {
00693        COLOR  ctmp;
00694        double d;
00695 
00696        d = xyYin[2] / xyYin[1];
00697        ctmp[0] = xyYin[0] * d;
00698        ctmp[1] = xyYin[2];
00699        ctmp[2] = (1. - xyYin[0] - xyYin[1]) * d;
00700                             /* allow negative values */
00701        colortrans(rgbout, xyz2rgbmat, ctmp);
00702 }
00703 
00704 
00705 static void
00706 picdebug(void)                     /* put out debugging picture */
00707 {
00708        static COLOR  blkcol = BLKCOLOR;
00709        COLOR  *scan;
00710        int    y, i;
00711        register int  x, rg;
00712 
00713        if (fseek(stdin, 0L, 0) == EOF) {
00714               fprintf(stderr, "%s: cannot seek on input picture\n", progname);
00715               exit(1);
00716        }
00717        getheader(stdin, NULL, NULL);             /* skip input header */
00718        fgetresolu(&xmax, &ymax, stdin);
00719                                           /* allocate scanline */
00720        scan = (COLOR *)malloc(xmax*sizeof(COLOR));
00721        if (scan == NULL) {
00722               perror(progname);
00723               exit(1);
00724        }
00725                                           /* finish debug header */
00726        fputformat(COLRFMT, debugfp);
00727        putc('\n', debugfp);
00728        fprtresolu(xmax, ymax, debugfp);
00729                                           /* write debug picture */
00730        for (y = ymax-1; y >= 0; y--) {
00731               if (freadscan(scan, xmax, stdin) < 0) {
00732                      fprintf(stderr, "%s: error rereading input picture\n",
00733                                    progname);
00734                      exit(1);
00735               }
00736               for (x = 0; x < xmax; x++) {
00737                      rg = chartndx(x, y, &i);
00738                      if (rg == RG_CENT) {
00739                             if (!(1L<<i & gmtflags) || (x+y)&07) {
00740                                    copycolor(scan[x], mbRGB[i]);
00741                                    clipgamut(scan[x], bright(scan[x]),
00742                                           CGAMUT, colmin, colmax);
00743                             } else
00744                                    copycolor(scan[x], blkcol);
00745                      } else if (rg == RG_CORR)
00746                             cvtcolor(scan[x], scan[x]);
00747                      else if (rg != RG_ORIG)
00748                             copycolor(scan[x], blkcol);
00749               }
00750               if (fwritescan(scan, xmax, debugfp) < 0) {
00751                      fprintf(stderr, "%s: error writing debugging picture\n",
00752                                    progname);
00753                      exit(1);
00754               }
00755        }
00756                                           /* clean up */
00757        fclose(debugfp);
00758        free((void *)scan);
00759 }
00760 
00761 
00762 static void
00763 clrdebug(void)                     /* put out debug picture from color input */
00764 {
00765        static COLR   blkclr = BLKCOLR;
00766        COLR   mbclr[24], cvclr[24], orclr[24];
00767        COLR   *scan;
00768        COLOR  ctmp, ct2;
00769        int    y, i;
00770        register int  x, rg;
00771                                           /* convert colors */
00772        for (i = 0; i < 24; i++) {
00773               copycolor(ctmp, mbRGB[i]);
00774               clipgamut(ctmp, bright(ctmp), CGAMUT, cblack, cwhite);
00775               setcolr(mbclr[i], colval(ctmp,RED),
00776                             colval(ctmp,GRN), colval(ctmp,BLU));
00777               if (inpflags & 1L<<i) {
00778                      copycolor(ctmp, inpRGB[i]);
00779                      clipgamut(ctmp, bright(ctmp), CGAMUT, cblack, cwhite);
00780                      setcolr(orclr[i], colval(ctmp,RED),
00781                                    colval(ctmp,GRN), colval(ctmp,BLU));
00782                      if (rawmap)
00783                             copycolr(cvclr[i], mbclr[i]);
00784                      else {
00785                             bresp(ctmp, inpRGB[i]);
00786                             colortrans(ct2, solmat, ctmp);
00787                             clipgamut(ct2, bright(ct2), CGAMUT,
00788                                           cblack, cwhite);
00789                             setcolr(cvclr[i], colval(ct2,RED),
00790                                           colval(ct2,GRN),
00791                                           colval(ct2,BLU));
00792                      }
00793               }
00794        }
00795                                           /* allocate scanline */
00796        scan = (COLR *)malloc(xmax*sizeof(COLR));
00797        if (scan == NULL) {
00798               perror(progname);
00799               exit(1);
00800        }
00801                                           /* finish debug header */
00802        fputformat(COLRFMT, debugfp);
00803        putc('\n', debugfp);
00804        fprtresolu(xmax, ymax, debugfp);
00805                                           /* write debug picture */
00806        for (y = ymax-1; y >= 0; y--) {
00807               for (x = 0; x < xmax; x++) {
00808                      rg = chartndx(x, y, &i);
00809                      if (rg == RG_CENT) {
00810                             if (!(1L<<i & gmtflags) || (x+y)&07)
00811                                    copycolr(scan[x], mbclr[i]);
00812                             else
00813                                    copycolr(scan[x], blkclr);
00814                      } else if (rg == RG_BORD || !(1L<<i & inpflags))
00815                             copycolr(scan[x], blkclr);
00816                      else if (rg == RG_ORIG)
00817                             copycolr(scan[x], orclr[i]);
00818                      else /* rg == RG_CORR */
00819                             copycolr(scan[x], cvclr[i]);
00820               }
00821               if (fwritecolrs(scan, xmax, debugfp) < 0) {
00822                      fprintf(stderr, "%s: error writing debugging picture\n",
00823                                    progname);
00824                      exit(1);
00825               }
00826        }
00827                                           /* clean up */
00828        fclose(debugfp);
00829        free((void *)scan);
00830 }
00831 
00832 
00833 static void
00834 getpos(              /* get boundary position */
00835        char   *name,
00836        int    bnds[2],
00837        FILE   *fp
00838 )
00839 {
00840        char   buf[64];
00841 
00842        fprintf(stderr, "\tSelect corner: %s\n", name);
00843        if (fgets(buf, sizeof(buf), fp) == NULL ||
00844                      sscanf(buf, "%d %d", &bnds[0], &bnds[1]) != 2) {
00845               fprintf(stderr, "%s: read error from display process\n",
00846                             progname);
00847               exit(1);
00848        }
00849 }
00850 
00851 
00852 static void
00853 pickchartpos(        /* display picture and pick chart location */
00854        char   *pfn
00855 )
00856 {
00857        char   combuf[PATH_MAX];
00858        FILE   *pfp;
00859 
00860        sprintf(combuf, DISPCOM, pfn);
00861        if ((pfp = popen(combuf, "r")) == NULL) {
00862               perror(combuf);
00863               exit(1);
00864        }
00865        fputs("Use middle mouse button to select chart corners:\n", stderr);
00866        getpos("upper left (dark skin)", bounds[0], pfp);
00867        getpos("upper right (bluish green)", bounds[1], pfp);
00868        getpos("lower left (white)", bounds[2], pfp);
00869        getpos("lower right (black)", bounds[3], pfp);
00870        fputs("Got it -- quit display program.\n", stderr);
00871        pclose(pfp);
00872 }