ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/macbethcal.c
(Generate patch)

Comparing ray/src/px/macbethcal.c (file contents):
Revision 2.1 by greg, Wed Oct 11 10:40:13 1995 UTC vs.
Revision 2.3 by greg, Wed Oct 11 18:52:00 1995 UTC

# Line 51 | Line 51 | 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 < #define NMBMOD          3       /* Number of MacBeth moderate colors */
55 < short   mbmod[NMBMOD] = {1,3,2};
54 > #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   #define NMBSAT          6       /* Number of MacBeth saturated colors */
57   short   mbsat[NMBSAT] = {14,12,13,15,16,17};
58  
# Line 62 | Line 62 | double imgxfm[3][3];           /* coordinate transformation mat
62  
63   COLOR   picRGB[24];             /* picture colors */
64  
65 < double  bramp[NMBNEU][2];       /* brightness ramp */
65 > COLOR   bramp[NMBNEU][2];       /* brightness ramp (per primary) */
66   double  solmat[3][3];           /* color mapping matrix */
67  
68   FILE    *debugfp = NULL;
# Line 125 | Line 125 | char   **argv;
125          }
126          init();                         /* initialize */
127          getcolors();                    /* get picture colors */
128 +        compute();                      /* compute color mapping */
129          putmapping();                   /* put out color mapping */
130          putdebug();                     /* put out debug picture */
131          exit(0);
# Line 138 | Line 139 | userr:
139   init()                          /* initialize */
140   {
141          double  quad[4][2];
141 –        COLOR   ctmp;
142 –        double  d;
143 –        int     i;
142                                          /* make coordinate transformation */
143          quad[0][0] = bounds[0][0];
144          quad[0][1] = bounds[0][1];
# Line 155 | Line 153 | init()                         /* initialize */
153                  fprintf(stderr, "%s: bad chart boundaries\n", progname);
154                  exit(1);
155          }
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 –        }
156   }
157  
158  
# Line 238 | Line 228 | getcolors()                            /* load in picture colors */
228   }
229  
230  
231 < double
232 < bresp(x)                /* piecewise linear interpolation of brightness */
243 < double  x;
231 > bresp(y, x)             /* piecewise linear interpolation of primaries */
232 > COLOR   y, x;
233   {
234 <        register int    n = NMBNEU;
234 >        double  cv[3];
235 >        register int    i, n;
236  
237 <        while (n > 0 && x < bramp[--n][0])
238 <                ;
239 <        return( ((bramp[n+1][0] - x)*bramp[n][1] +
240 <                                (x - bramp[n][0])*bramp[n+1][1]) /
241 <                        (bramp[n+1][0] - bramp[n][0]) );
237 >        for (i = 0; i < 3; i++) {
238 >                n = NMBNEU;
239 >                while (n > 0 && colval(x,i) < colval(bramp[--n][0],i))
240 >                        ;
241 >                cv[i] = ((colval(bramp[n+1][0],i) - colval(x,i)) *
242 >                                                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 >                if (cv[i] < 0.) cv[i] = 0.;
247 >        }
248 >        setcolor(y, cv[0], cv[1], cv[2]);
249   }
250  
251  
252 < putmapping()                    /* compute and print mapping for pcomb -f */
252 > compute()                       /* compute color mapping */
253   {
254 <        float   clrin[NMBMOD][3], clrout[NMBMOD][3];
255 <        register int    i, j;
254 >        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                                          /* compute piecewise luminance curve */
267          for (i = 0; i < NMBNEU; i++) {
268 <                bramp[i][0] = bright(picRGB[mbneu[i]]);
269 <                bramp[i][1] = bright(mbRGB[mbneu[i]]);
268 >                copycolor(bramp[i][0], picRGB[mbneu[i]]);
269 >                copycolor(bramp[i][1], mbRGB[mbneu[i]]);
270          }
271                                          /* compute color matrix */
272 <        for (i = 0; i < NMBMOD; i++)
273 <                for (j = 0; j < 3; j++) {
274 <                        clrin[i][j] = bresp(picRGB[mbmod[i]][j]);
275 <                        clrout[i][j] = mbRGB[mbmod[i]][j];
269 <                }
272 >        for (i = 0; i < NMBMOD; i++) {
273 >                bresp(clrin[i], picRGB[mbmod[i]]);
274 >                copycolor(clrout[i], mbRGB[mbmod[i]]);
275 >        }
276          compsoln(clrin, clrout, NMBMOD);
277 + }
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                                          /* print brightness mapping */
285 <        printf("xval(i) : select(i, 0");
286 <        for (i = 0; i < NMBNEU; i++)
287 <                printf(", %g", bramp[i][0]);
288 <        printf(");\n");
289 <        printf("yval(i) : select(i, 0");
290 <        for (i = 0; i < NMBNEU; i++)
291 <                printf(", %g", bramp[i][1]);
292 <        printf(");\n");
293 <        printf("ifind(x,f,n) : if(1.5-n, 1, if(x-f(n), n, ifind(x,f,n-1)));\n");
294 <        printf("binterp(i,x) : ((xval(i+1)-x)*yval(i) + (x-xval(i))*yval(i+1))/\n");
295 <        printf("\t\t(xval(i+1) - xval(i));\n");
296 <        printf("bresp(x) : binterp(ifind(x,xval,xval(0)-1), x);\n");
297 <        printf("nred = bresp(ri(1));\n");
298 <        printf("ngrn = bresp(gi(1));\n");
299 <        printf("nblu = bresp(bi(1));\n");
285 >        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                                          /* print color mapping */
309 <        printf("ro = %g*nred + %g*ngrn + %g*nblu\n",
310 <                        solmat[0][0], solmat[1][0], solmat[2][0]);
311 <        printf("go = %g*nred + %g*ngrn + %g*nblu\n",
312 <                        solmat[0][1], solmat[1][1], solmat[2][1]);
313 <        printf("bo = %g*nred + %g*ngrn + %g*nblu\n",
314 <                        solmat[0][2], solmat[1][2], solmat[2][2]);
309 >        printf("ro = %g*rn + %g*gn + %g*bn ;\n",
310 >                        solmat[0][0], solmat[0][1], solmat[0][2]);
311 >        printf("go = %g*rn + %g*gn + %g*bn ;\n",
312 >                        solmat[1][0], solmat[1][1], solmat[1][2]);
313 >        printf("bo = %g*rn + %g*gn + %g*bn ;\n",
314 >                        solmat[2][0], solmat[2][1], solmat[2][2]);
315   }
316  
317  
318 + #if NMBMOD == 3
319   compsoln(cin, cout, n)          /* solve 3x3 system */
320 < float   cin[][3], cout[][3];
320 > COLOR   cin[], cout[];
321   int     n;
322   {
323          extern double   mx3d_adjoint(), fabs();
# Line 310 | Line 332 | int    n;
332          }
333          for (i = 0; i < 3; i++)
334                  for (j = 0; j < 3; j++)
335 <                        mat[i][j] = cin[j][i];
335 >                        mat[i][j] = colval(cin[j],i);
336          det = mx3d_adjoint(mat, invmat);
337          if (fabs(det) < 1e-4) {
338                  fprintf(stderr, "%s: cannot compute color mapping\n",
# Line 325 | Line 347 | int    n;
347                          invmat[i][j] /= det;
348          for (i = 0; i < 3; i++) {
349                  for (j = 0; j < 3; j++)
350 <                        rowv[j] = cout[j][i];
351 <                mx3d_transform(rowv, invmat, colv);
350 >                        colv[j] = colval(cout[j],i);
351 >                mx3d_transform(colv, invmat, rowv);
352                  for (j = 0; j < 3; j++)
353 <                        solmat[j][i] = colv[j];
353 >                        solmat[i][j] = rowv[j];
354          }
355   }
356 + #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  
365 +        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  
402 +
403   cvtcolor(cout, cin)             /* convert color according to our mapping */
404   COLOR   cout, cin;
405   {
406          double  r, g, b;
340 –        double  r1, g1, b1;
407  
408 <        r = bresp(colval(cin,RED));
409 <        g = bresp(colval(cin,GRN));
410 <        b = bresp(colval(cin,BLU));
411 <        r1 = r*solmat[0][0] + g*solmat[1][0] + b*solmat[2][0];
412 <        if (r1 < 0) r1 = 0;
413 <        g1 = r*solmat[0][1] + g*solmat[1][1] + b*solmat[2][1];
414 <        if (g1 < 0) g1 = 0;
415 <        b1 = r*solmat[0][2] + g*solmat[1][2] + b*solmat[2][2];
416 <        if (b1 < 0) b1 = 0;
417 <        setcolor(cout, r1, g1, b1);
408 >        bresp(cout, cin);
409 >        r = colval(cout,0)*solmat[0][0] + colval(cout,1)*solmat[0][1]
410 >                        + colval(cout,2)*solmat[0][2];
411 >        if (r < 0) r = 0;
412 >        g = colval(cout,0)*solmat[1][0] + colval(cout,1)*solmat[1][1]
413 >                        + colval(cout,2)*solmat[1][2];
414 >        if (g < 0) g = 0;
415 >        b = colval(cout,0)*solmat[2][0] + colval(cout,1)*solmat[2][1]
416 >                        + colval(cout,2)*solmat[2][2];
417 >        if (b < 0) b = 0;
418 >        setcolor(cout, r, g, b);
419   }
420  
421  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines