home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / cprog / illum.lha / clr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-01-03  |  29.6 KB  |  883 lines
  1. /* ****************************************************************
  2.  *                           clr.c
  3.  * ****************************************************************
  4.  *  MODULE PURPOSE:
  5.  *      This module contains routines for color
  6.  *      transformations, color space manipulations,
  7.  *      and spectral sampling.
  8.  *
  9.  *  MODULE CONTENTS:
  10.  *      CLR_init            - initialized the color routines
  11.  *      CLR_read_mtl        - read a material file
  12.  *      CLR_blackbody       - returns the spectral curve for
  13.  *                              a blackbody at the specified
  14.  *                              temperature.
  15.  *      CLR_add_spect       - add two spectral curves
  16.  *      CLR_mult_spect      - multiply two spectral curves
  17.  *      CLR_scale_spect     - scale a spectral curve
  18.  *      CLR_area_spect      - get the area under a curve
  19.  *      CLR_spect_to_xyz    - sample a curve to xyz
  20.  *      CLR_spect_to_rgb    - sample a curve to rgb
  21.  *      CLR_get_CIEXYZ      - get the CIEXYZ matching functions
  22.  *      CLR_get_rgb         - returns the color space primaries
  23.  *      CLR_get_min_wl      - returns the max wavelength
  24.  *      CLR_get_max_wl      - returns the min wavelength
  25.  *      CLR_get_xyz_rgb     - returns the xyz to rgb matrix
  26.  *      CLR_get_rgb_xyz     - returns the rgb to xyz matrix
  27.  *      CLR_get_rgb_yiq     - returns the rgb to yiq matrix
  28.  *      CLR_get_yiq_rgb     - returns the yiq to rgb matrix
  29.  *      CLR_rgb_to_aux_rgb  - returns the matrix from the
  30.  *                              current color space to an
  31.  *                              auxiliary color space
  32.  *      CLR_xyz_lab         - transforms from xyz to Lab
  33.  *      CLR_xyz_luv         - transforms from xyz to Luv
  34.  *      CLR_t_concat        - concatenate color transforms
  35.  *      CLR_t_inverse       - return an inverse transform
  36.  *      CLR_exit            - finish with the color routines
  37.  *
  38.  *  NOTES:
  39.  *      > Spectral curve manipulations use the lower and upper
  40.  *          bounds established at init.  Curves are at 1nm
  41.  *          increments and are arrays of doubles, (max - min + 1)
  42.  *          elements long.
  43.  *      > Matrix scaling is so that an RGB of 1,1,1 transforms
  44.  *          to an XYZ with a Y value of 1.0
  45.  *      > The XYZ and RGB sampling from spectral curves are
  46.  *          scaled so that that a mirror reflector (1.0 for all
  47.  *          wavelengths) has a Y value of 1.0.
  48.  *      > Wavelength bounds of 380nm to 780nm (the visible
  49.  *          spectrum) are appropriate for most applications.
  50.  */
  51. #include <stdio.h>
  52. #include <math.h>
  53. #include "clr.h"
  54.  
  55. #define RED         0
  56. #define GREEN       1
  57. #define BLUE        2
  58. #define WHITE       3
  59. #define LOAD_MAT(a,b) \
  60. {   int ii,jj; \
  61.     for (ii=0; ii<=2; ii++) \
  62.     for (jj=0; jj<=2; jj++) a[ii][jj] = b[ii][jj]; \
  63. }
  64.  
  65. static int          init = FALSE;
  66. static int          min_wl = 380;
  67. static int          max_wl = 780;
  68. static double       xyz_rgb_mat[3][3];
  69. static double       rgb_xyz_mat[3][3];
  70. static double       *X_tristim = NULL;
  71. static double       *Y_tristim = NULL;
  72. static double       *Z_tristim = NULL;
  73. static double       *work_curve = NULL;
  74. static double       xyz_scale = 1.0;
  75. static CLR_XYZ      rgb_primary[4];
  76. static int          clr__cspace_to_xyz ();
  77. static CLR_XYZ      white;
  78. static CLR_LUV      ref_luv;
  79.  
  80. /* This is the RGB to YIQ matrix and YIQ to RGB matrix as
  81.  *  given in Sect 5.2 NTSC and RGB Video
  82.  */
  83. static double       rgb_yiq_mat[3][3] = {
  84.     {0.299,     0.587,      0.144},
  85.     {0.596,    -0.275,     -0.321},
  86.     {0.212,    -0.528,      0.311}};
  87.  
  88. static double       yiq_rgb_mat[3][3] = {
  89.     {1.0000,    0.9557,     0.6199},
  90.     {1.0000,   -0.2716,    -0.6469},
  91.     {1.0000,   -1.1082,     1.7051}};
  92.  
  93. /* This is the NTSC primaries with D6500 white point for use as
  94.  *  the default initialization as given in sect 5.1.1 Color
  95.  *  Correction for Display.
  96.  */
  97. static CLR_XYZ      rgb_NTSC[4] = {
  98.     {0.670,     0.330,      0.000},     /* red */
  99.     {0.210,     0.710,      0.080},     /* green */
  100.     {0.140,     0.080,      0.780},     /* blue */
  101.     {0.313,     0.329,      0.358}};    /* white */
  102.  
  103. /* This is the data for the CIEXYZ curves take from Judd and
  104.  *  Wyszecki (1975), table 2.6, these are for the 1931 standard
  105.  *  observer with a 2-degree visual field.
  106.  */
  107. static double       CIEXYZ[81][4] = {
  108.     {380, 0.0014, 0.0000, 0.0065}, {385, 0.0022, 0.0001, 0.0105},
  109.     {390, 0.0042, 0.0001, 0.0201}, {395, 0.0076, 0.0002, 0.0362},
  110.     {400, 0.0143, 0.0004, 0.0679}, {405, 0.0232, 0.0006, 0.1102},
  111.     {410, 0.0435, 0.0012, 0.2074}, {415, 0.0776, 0.0022, 0.3713},
  112.     {420, 0.1344, 0.0040, 0.6456}, {425, 0.2148, 0.0073, 1.0391},
  113.     {430, 0.2839, 0.0116, 1.3856}, {435, 0.3285, 0.0168, 1.6230},
  114.     {440, 0.3483, 0.0230, 1.7471}, {445, 0.3481, 0.0298, 1.7826},
  115.     {450, 0.3362, 0.0380, 1.7721}, {455, 0.3187, 0.0480, 1.7441},
  116.     {460, 0.2908, 0.0600, 1.6692}, {465, 0.2511, 0.0739, 1.5281},
  117.     {470, 0.1954, 0.0910, 1.2876}, {475, 0.1421, 0.1126, 1.0419},
  118.     {480, 0.0956, 0.1390, 0.8310}, {485, 0.0580, 0.1693, 0.6162},
  119.     {490, 0.0320, 0.2080, 0.4652}, {495, 0.0147, 0.2586, 0.3533},
  120.     {500, 0.0049, 0.3230, 0.2720}, {505, 0.0024, 0.4073, 0.2123},
  121.     {510, 0.0093, 0.5030, 0.1582}, {515, 0.0291, 0.6082, 0.1117},
  122.     {520, 0.0633, 0.7100, 0.0782}, {525, 0.1096, 0.7932, 0.0573},
  123.     {530, 0.1655, 0.8620, 0.0422}, {535, 0.2257, 0.9149, 0.0298},
  124.     {540, 0.2904, 0.9540, 0.0203}, {545, 0.3597, 0.9803, 0.0134},
  125.     {550, 0.4334, 0.9950, 0.0087}, {555, 0.5121, 1.0000, 0.0057},
  126.     {560, 0.5945, 0.9950, 0.0039}, {565, 0.6784, 0.9786, 0.0027},
  127.     {570, 0.7621, 0.9520, 0.0021}, {575, 0.8425, 0.9154, 0.0018},
  128.     {580, 0.9163, 0.8700, 0.0017}, {585, 0.9786, 0.8163, 0.0014},
  129.     {590, 1.0263, 0.7570, 0.0011}, {595, 1.0567, 0.6949, 0.0010},
  130.     {600, 1.0622, 0.6310, 0.0008}, {605, 1.0456, 0.5668, 0.0006},
  131.     {610, 1.0026, 0.5030, 0.0003}, {615, 0.9384, 0.4412, 0.0002},
  132.     {620, 0.8544, 0.3810, 0.0002}, {625, 0.7514, 0.3210, 0.0001},
  133.     {630, 0.6424, 0.2650, 0.0000}, {635, 0.5419, 0.2170, 0.0000},
  134.     {640, 0.4479, 0.1750, 0.0000}, {645, 0.3608, 0.1382, 0.0000},
  135.     {650, 0.2835, 0.1070, 0.0000}, {655, 0.2187, 0.0816, 0.0000},
  136.     {660, 0.1649, 0.0610, 0.0000}, {665, 0.1212, 0.0446, 0.0000},
  137.     {670, 0.0874, 0.0320, 0.0000}, {675, 0.0636, 0.0232, 0.0000},
  138.     {680, 0.0468, 0.0170, 0.0000}, {685, 0.0329, 0.0119, 0.0000},
  139.     {690, 0.0227, 0.0082, 0.0000}, {695, 0.0158, 0.0057, 0.0000},
  140.     {700, 0.0114, 0.0041, 0.0000}, {705, 0.0081, 0.0029, 0.0000},
  141.     {710, 0.0058, 0.0021, 0.0000}, {715, 0.0041, 0.0015, 0.0000},
  142.     {720, 0.0029, 0.0010, 0.0000}, {725, 0.0020, 0.0007, 0.0000},
  143.     {730, 0.0014, 0.0005, 0.0000}, {735, 0.0010, 0.0004, 0.0000},
  144.     {740, 0.0007, 0.0002, 0.0000}, {745, 0.0005, 0.0002, 0.0000},
  145.     {750, 0.0003, 0.0001, 0.0000}, {755, 0.0002, 0.0001, 0.0000},
  146.     {760, 0.0002, 0.0001, 0.0000}, {765, 0.0001, 0.0000, 0.0000},
  147.     {770, 0.0001, 0.0000, 0.0000}, {775, 0.0001, 0.0000, 0.0000},
  148.     {780, 0.0000, 0.0000, 0.0000} };
  149.  
  150.  
  151.  
  152. /* ****************************************************************
  153.  * CLR_init (min, max, rgb)
  154.  *  int         min, max    (in) wavelength bounds (nm)
  155.  *  CLR_XYZ     rgb[4]      (in) the primaries.  If NULL,
  156.  *                              then the NTSC primaries with
  157.  *                              a D6500 white point are used.
  158.  *
  159.  * Initializes the color routines for a spectral range from
  160.  *  'min'nm to 'max'nm and an RGB color space given by 'rgb'
  161.  *  Returns TRUE if successful and FALSE on error.
  162.  */
  163. int CLR_init (min, max, rgb)
  164. int         min;
  165. int         max;
  166. CLR_XYZ     rgb[4];
  167. {   int     color, ct;
  168.     double  len;
  169.     char    *malloc();
  170.  
  171.     if (init) return FALSE;
  172.     init = TRUE;
  173.     min_wl = min;
  174.     max_wl = max;
  175.  
  176.     /* load primaries and build transformations, use the
  177.      * defaults is rgb==NULL
  178.      */
  179.     for (color=RED; color<=WHITE; color++) {
  180.     if (rgb == NULL) rgb_primary[color] = rgb_NTSC[color];
  181.     else rgb_primary[color] = rgb[color];
  182.     rgb_primary[color].z = 1.0 -
  183.         rgb_primary[color].x - rgb_primary[color].y;
  184.     }
  185.     if (!clr__cspace_to_xyz(rgb_primary, rgb_xyz_mat)) goto error;
  186.     if (!CLR_t_inverse(rgb_xyz_mat, xyz_rgb_mat)) goto error;
  187.  
  188.     /* build reference info for L*a*b*, L*u*v* transforms */
  189.     white.x =
  190.     (100.0 * rgb_primary[WHITE].x) / rgb_primary[WHITE].y;
  191.     white.y = 100.0;
  192.     white.z =
  193.     (100.0 * rgb_primary[WHITE].z) / rgb_primary[WHITE].y;
  194.     ref_luv.u = (4.0 * white.x) /
  195.     (white.x + (15.0 * white.y) + (3.0 * white.z));
  196.     ref_luv.v = (9.0 * white.y) /
  197.     (white.x + (15.0 * white.y) + (3.0 * white.z));
  198.  
  199.     /* build the XYZ sampling curves - allocate space for the
  200.      *  curves and interpolate the CIEXYZ table into continuous
  201.      *  curves at 1nm increments.
  202.      */
  203.     {   int     ii, nm;
  204.     double  *x, *y, *z;
  205.     double  x_cur, y_cur, z_cur;
  206.     double  x_inc, y_inc, z_inc;
  207.     if ((work_curve = (double *)malloc((unsigned)
  208.         (sizeof(double) * (max_wl-min_wl+1))))
  209.         == NULL) goto error;
  210.     if ((x = X_tristim = (double *)malloc((unsigned)
  211.         (sizeof(double) * (max_wl-min_wl+1))))
  212.         == NULL) goto error;
  213.     if ((y = Y_tristim = (double *)malloc((unsigned)
  214.         (sizeof(double) * (max_wl-min_wl+1))))
  215.         == NULL) goto error;
  216.     if ((z = Z_tristim = (double *)malloc((unsigned)
  217.         (sizeof(double) * (max_wl-min_wl+1))))
  218.         == NULL) goto error;
  219.     for (nm=min_wl; nm<380; nm++) *x++ = *y++ = *z++ = 0.0;
  220.     for (ct=0; ct<80; ct++) {
  221.         x_cur = CIEXYZ[ct][1];
  222.         y_cur = CIEXYZ[ct][2];
  223.         z_cur = CIEXYZ[ct][3];
  224.         x_inc = (CIEXYZ[ct+1][1] - x_cur) / 5.0;
  225.         y_inc = (CIEXYZ[ct+1][2] - y_cur) / 5.0;
  226.         z_inc = (CIEXYZ[ct+1][3] - z_cur) / 5.0;
  227.         for (ii=0; ii<5; ii++, nm++) {
  228.         if (nm > max_wl) goto XYZ_done;
  229.         if (nm >= min_wl) {
  230.             *x++ = x_cur; *y++ = y_cur; *z++ = z_cur;
  231.         }
  232.         x_cur += x_inc; y_cur += y_inc; z_cur += z_inc;
  233.         }
  234.     }
  235.     if (nm <= max_wl) {
  236.         *x++ = x_cur; *y++ = y_cur; *z++ = z_cur; nm++;
  237.     }
  238.     for ( ;nm<=max_wl; nm++) *x++ = *y++ = *z++ = 0.0;
  239.     }
  240. XYZ_done:
  241.     /* determine the scaling factor to be used in sampling spectral
  242.      *  curves to bring Y os a sampled identity curve to 1.
  243.      */
  244.     xyz_scale = 1.0 / CLR_area_spect (Y_tristim);
  245.  
  246.     return TRUE;
  247.  
  248. error:
  249.     CLR_exit();
  250.     return FALSE;
  251. }
  252.  
  253. /* ****************************************************************
  254.  * CLR_read_mtl (file, conduct, n, k, curve)
  255.  *  char        *file       (in) material file name
  256.  *  int         *conduct    (out) conductor - set to FALSE if not
  257.  *                              specified in the material file
  258.  *  double      *n, *k      (out) n and k, - set to 0.0 if not
  259.  *                              specified in the material file
  260.  *  double      *curve      (out) curve array
  261.  *
  262.  * Returns TRUE if the material was read, FALSE if the material
  263.  *  could not be found or an error was detected in the material
  264.  *  file.  NULL pointers can be given as arguments for 'conduct',
  265.  *  'n', 'k', and/or 'curve' if the properties are not of interest.
  266.  *
  267.  * The material file must be in the format given in Sect III.7.1
  268.  *  Material Data.
  269.  */
  270. int CLR_read_mtl (file, conduct, n, k, curve)
  271. char        *file;
  272. int         *conduct;
  273. double      *n, *k;
  274. double      *curve;
  275. {   FILE        *fp, *fopen();
  276.     char        in_str[200];
  277.     int         flag = 1;
  278.     int         cur_wl, last_wl, ind_wl;
  279.     double      cur_val, last_val, inc_val;
  280.  
  281.     if (!init) return FALSE;
  282.     if ((fp=fopen(file,"r")) == NULL) return FALSE;
  283.     if (n != NULL) *n = 0.0;
  284.     if (k != NULL) *k = 0.0;
  285.     if (conduct != NULL) *conduct = FALSE;
  286.     if (curve == NULL) flag = -1;
  287.     ind_wl = min_wl;
  288.     while (fgets(in_str,200,fp)) {
  289.     switch (in_str[0]) {
  290.         case '#':       /* comment */
  291.         break;
  292.         case 'k':       /* absorption coefficient */
  293.         if (k != NULL) (void)sscanf (in_str, "k %lf", k);
  294.         break;
  295.         case 'n':       /* index of refraction */
  296.         if (n != NULL) (void)sscanf (in_str, "n %lf", n);
  297.         break;
  298.         case 'c':       /* conductor */
  299.         if (conduct != NULL) *conduct = TRUE;
  300.         break;
  301.         case 'd':       /* dielectric */
  302.         if (conduct != NULL) *conduct = FALSE;
  303.         break;
  304.         default:        /* spectral data */
  305.         if ((flag != -1) && (sscanf(in_str,
  306.             "%d %lf", &cur_wl, &cur_val) == 2)) {
  307.             if (flag == 1) {
  308.             flag = 0;
  309.             while (ind_wl <= cur_wl) {
  310.                 *curve++ = cur_val;
  311.                 if (++ind_wl > max_wl) {
  312.                 flag == -1;
  313.                 break;
  314.                 }
  315.             }
  316.             }
  317.             else {
  318.             if (cur_wl < last_wl) {
  319.                 (void)fclose(fp);
  320.                 return FALSE;
  321.             }
  322.             inc_val = (cur_val - last_val) /
  323.                 (cur_wl - last_wl);
  324.             while (last_wl < ind_wl) {
  325.                 last_val += inc_val;
  326.                 last_wl++;
  327.             }
  328.             while (ind_wl <= cur_wl) {
  329.                 *curve++ = last_val;
  330.                 last_val += inc_val;
  331.                 if (++ind_wl > max_wl) {
  332.                 flag = -1;
  333.                 break;
  334.                 }
  335.             }
  336.             }
  337.             last_wl = cur_wl;
  338.             last_val = cur_val;
  339.         }
  340.     }
  341.     }
  342.     while (ind_wl <= max_wl) {
  343.     *curve++ = last_val;
  344.     ind_wl++;
  345.     }
  346.     (void)fclose(fp);
  347.     return TRUE;
  348. }
  349.  
  350. /* ****************************************************************
  351.  * CLR_blackbody (temp, curve)
  352.  *  double      temp        (in) temperature (degrees K)
  353.  *  double      *curve      (out) curve array
  354.  *
  355.  * Fills the 'curve' array with the illumination values for an
  356.  *  ideal blackbody radiator at the specified temperature.
  357.  *  Following the convention used for normalizing other illuminant
  358.  *  curves, the blackbody curve is normalized so that the value
  359.  *  at 560nm is 1.0.  The formulation is taken from Judd and
  360.  *  Wysecki (1975) p.106.  NOTE: c1 and c2 have been adjusted
  361.  *  to account for the wavelength to be in nm.
  362.  */
  363. int CLR_blackbody (temp, curve)
  364. double  temp;
  365. double  *curve;
  366. {   double  norm;
  367.     double  c1 = 3.74150E29;
  368.     double  c2 = 1.4388E7;
  369.     double  wl;
  370.     int     ct;
  371.  
  372.     if (temp < 0.0) return FALSE;
  373.     norm = (pow((double)560.0, (double)5.0) *
  374.     (exp (c2 / (560.0 * temp)) - 1.0)) / c1;
  375.     for (ct=max_wl-min_wl+1, wl=min_wl; --ct>=0; wl += 1.0) {
  376.     *curve++ = (norm * c1) / (pow(wl, (double)5.0) *
  377.         (exp (c2 / (wl * temp)) - 1.0));
  378.     }
  379.     return TRUE;
  380. }
  381.  
  382. /* ****************************************************************
  383.  * CLR_add_spect (c1, c2, c3)
  384.  *  double  *c1, *c2    (in)  - input curves
  385.  *  double  *c3         (out) - output curve
  386.  *
  387.  * Add spectral curve 'c1' to spectral curve 'c2' to get spectral
  388.  *  curve 'c3'.  Returns TRUE if successful, FALSE if the CLR_
  389.  *  routines are not initialized.
  390.  */
  391. CLR_add_spect (c1,c2,c3)
  392. double  *c1, *c2, *c3;
  393. {   int     ct;
  394.     if (!init) return FALSE;
  395.     for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ + *c2++;
  396.     return TRUE;
  397. }
  398.  
  399. /* ****************************************************************
  400.  * CLR_mult_spect (c1, c2, c3)
  401.  *  double  *c1, *c2    (in)  - input curves
  402.  *  double  *c3         (out) - output curve
  403.  *
  404.  * Multiply spectral curve 'c1' by spectral curve 'c2' to get
  405.  *  spectral curve 'c3'.  Returns TRUE if successful, FALSE if
  406.  *  the CLR_ routines are not initialized.
  407.  */
  408. CLR_mult_spect (c1,c2,c3)
  409. double  *c1, *c2, *c3;
  410. {   int     ct;
  411.     if (!init) return FALSE;
  412.     for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ * *c2++;
  413.     return TRUE;
  414. }
  415.  
  416. /* ****************************************************************
  417.  * CLR_scale_spect (c1, scale, c3)
  418.  *  double  *c1         (in)  - input curve
  419.  *  double  scale       (in)  - scale
  420.  *  double  *c3         (out) - output curve
  421.  *
  422.  * Multiply spectral curve 'c1' by 'scale' to get spectral curve
  423.  *  'c3'.  Returns TRUE if successful, FALSE if the CLR_ routines
  424.  *  are not initialized.
  425.  */
  426. CLR_scale_spect (c1,scale,c3)
  427. double  *c1, scale, *c3;
  428. {   int     ct;
  429.     if (!init) return FALSE;
  430.     for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ * scale;
  431.     return TRUE;
  432. }
  433.  
  434. /* ****************************************************************
  435.  * double CLR_area_spect (c1)
  436.  *  double  *c1         (in)  - input curve
  437.  *
  438.  * Returns the are under the spectral curve 'c1'.  Returns 0 if
  439.  *  the CLR_ routines are not initialized.
  440.  */
  441. double CLR_area_spect (c1)
  442. double  *c1;
  443. {   int     ct;
  444.     double  area = 0.0;
  445.     if (!init) return FALSE;
  446.     for (ct=max_wl-min_wl+1; --ct>=0; ) area += *c1++;
  447.     return area;
  448. }
  449.  
  450. /* ****************************************************************
  451.  * CLR_XYZ CLR_spect_to_xyz(spectral)
  452.  *  double  *spectral   (in) - the spectral curve
  453.  *
  454.  * Returns the sample values if successful, and a sample value
  455.  *  of (0,0,0) if the CLR_ routines are not initialized.
  456.  *
  457.  * Multiplies the spectral curve by each of the sampling curves
  458.  *  then integrates the resulting curves.  The XYZ values are then
  459.  *  normalized such that an identity material has Y=1.
  460.  */
  461. CLR_XYZ CLR_spect_to_xyz (spectral)
  462. double      *spectral;
  463. {   CLR_XYZ     xyz;
  464.     if (!init) {xyz.x = xyz.y = xyz.z = 0.0; return xyz;}
  465.     (void)CLR_mult_spect(spectral, X_tristim, work_curve);
  466.     xyz.x = xyz_scale * CLR_area_spect(work_curve);
  467.     (void)CLR_mult_spect(spectral, Y_tristim, work_curve);
  468.     xyz.y = xyz_scale * CLR_area_spect(work_curve);
  469.     (void)CLR_mult_spect(spectral, Z_tristim, work_curve);
  470.     xyz.z = xyz_scale * CLR_area_spect(work_curve);
  471.     return xyz;
  472. }
  473.  
  474. /* ****************************************************************
  475.  * CLR_RGB CLR_spect_to_rgb (spectral)
  476.  *  double  *spectral   (in) - the spectral curve
  477.  *
  478.  * Returns the sample values if successful, and a sample value of
  479.  *  (0,0,0) if the CLR_ routines are not initialized.
  480.  *
  481.  * The curve is first sampled to XYZ then transformed to RGB
  482.  */
  483. CLR_RGB CLR_spect_to_rgb (spectral)
  484. double      *spectral;
  485. {   CLR_XYZ     xyz;
  486.     CLR_RGB     rgb;
  487.     if (!init) {rgb.r = rgb.g = rgb.b = 0.0; return rgb;}
  488.     xyz = CLR_spect_to_xyz (spectral);
  489.     rgb.r = (xyz_rgb_mat[0][0] * xyz.x) +
  490.     (xyz_rgb_mat[0][1] * xyz.y) + (xyz_rgb_mat[0][2] * xyz.z);
  491.     rgb.g = (xyz_rgb_mat[1][0] * xyz.x) +
  492.     (xyz_rgb_mat[1][1] * xyz.y) + (xyz_rgb_mat[1][2] * xyz.z);
  493.     rgb.b = (xyz_rgb_mat[2][0] * xyz.x) +
  494.     (xyz_rgb_mat[2][1] * xyz.y) + (xyz_rgb_mat[2][2] * xyz.z);
  495.     return rgb;
  496. }
  497.  
  498. /* ****************************************************************
  499.  * CLR_get_CIEXYZ (X,Y,Z)
  500.  *  double  *X,*Y,*Z    (mod) arrays to hold the curves
  501.  *
  502.  * Copies the XYZ sampling curves into the user supplied buffers.
  503.  *  Returns TRUE is successful and FALSE if the CLR_ routines are
  504.  *  not initialized.
  505.  */
  506. CLR_get_CIEXYZ(X,Y,Z)
  507. double      *X,*Y,*Z;
  508. {   int     ct;
  509.     double  *x,*y,*z;
  510.     if (!init) return FALSE;
  511.     x = X_tristim; y = Y_tristim; z = Z_tristim;
  512.     for (ct=max_wl-min_wl+1; --ct>=0; ) {
  513.     *X++ = *x++; *Y++ = *y++; *Z++ = *z++;
  514.     }
  515.     return TRUE;
  516. }
  517.  
  518. /* ****************************************************************
  519.  * CLR_get_rgb (rgb)
  520.  *  CLR_XYZ     rgb[4]  (mod) the primaries
  521.  *
  522.  * Fills 'rgb' with the primaries the CLR_ routines are initialized
  523.  *  to.  Returns TRUE if successful and FALSE if the CLR_ routines
  524.  *  are not initialized.
  525.  */
  526. CLR_get_rgb (rgb)
  527. CLR_XYZ     *rgb;
  528. {   int     ct;
  529.     if (!init) return FALSE;
  530.     for (ct=0; ct<=3; ct++) rgb[ct] = rgb_primary[ct];
  531.     return TRUE;
  532. }
  533.  
  534. /* ****************************************************************
  535.  * CLR_get_min_wl ()
  536.  *
  537.  * Returns the minimum wavelength bound for which the CLR_ routines
  538.  *  are initialized, returns -1 if the CLR_ routines are not
  539.  *  initialized.
  540.  */
  541. CLR_get_min_wl()
  542. {   if (!init) return -1;
  543.     return min_wl;
  544. }
  545.  
  546. /* ****************************************************************
  547.  * CLR_get_max_wl ()
  548.  *
  549.  * Returns the minimum wavelength bound for which the CLR_ routines
  550.  *  are initialized, returns -1 if the CLR_ routines are not
  551.  *  initialized.
  552.  */
  553. CLR_get_max_wl()
  554. {   if (!init) return -1;
  555.     return max_wl;
  556. }
  557.  
  558. /* ****************************************************************
  559.  * CLR_get_xyz_rgb (mat)
  560.  *  double        mat[3][3]     (mod) - matrix to be loaded
  561.  *
  562.  *  Copies the CIEXYZ to RGB transformation into the user supplied
  563.  *   matrix.  Returns TRUE if successful and FALSE if the CLR_
  564.  *   routines are not initialized.
  565.  */
  566. int CLR_get_xyz_rgb(mat)
  567. double mat[][3];
  568. {   if (!init) return FALSE;
  569.     LOAD_MAT(mat, xyz_rgb_mat);
  570.     return TRUE;
  571. }
  572.  
  573. /* ****************************************************************
  574.  * CLR_get_rgb_xyz (mat)
  575.  *  double        mat[3][3]     (mod) - matrix to be loaded
  576.  *
  577.  *  Copies the RGB to CIEXYZ transformation into the user supplied
  578.  *   matrix.  Returns TRUE if successful and FALSE if the CLR_
  579.  *   routines are not initialized.
  580.  */
  581. int CLR_get_rgb_xyz(mat)
  582. double  mat[][3];
  583. {   if (!init) return FALSE;
  584.     LOAD_MAT(mat, rgb_xyz_mat);
  585.     return TRUE;
  586. }
  587.  
  588. /* ****************************************************************
  589.  * CLR_get_yiq_rgb (mat)
  590.  *  double        mat[3][3]     (mod) - matrix to be loaded
  591.  *
  592.  *  Copies the YIQ to RGB transformation into the user supplied
  593.  *   matrix.  Returns TRUE if successful and FALSE if the CLR_
  594.  *   routines are not initialized.
  595.  */
  596. int CLR_get_yiq_rgb(mat)
  597. double mat[][3];
  598. {   if (!init) return FALSE;
  599.     LOAD_MAT(mat, yiq_rgb_mat);
  600.     return TRUE;
  601. }
  602.  
  603. /* ****************************************************************
  604.  * CLR_get_rgb_yiq (mat)
  605.  *  double        mat[3][3]     (mod) - matrix to be loaded
  606.  *
  607.  *  Copies the RGB to YIQ transformation into the user supplied
  608.  *   matrix.  Returns TRUE if successful and FALSE if the CLR_
  609.  *   routines are not initialized.
  610.  */
  611. int CLR_get_rgb_yiq(mat)
  612. double  mat[][3];
  613. {   if (!init) return FALSE;
  614.     LOAD_MAT(mat, rgb_yiq_mat);
  615.     return TRUE;
  616. }
  617.  
  618. /* ****************************************************************
  619.  * CLR_rgb_to_aux_rgb (to, from, rgb_aux)
  620.  *  double        to[3][3]      (mod) - matrix to aux rgb
  621.  *  double        from[3][3]    (mod) - matrix from aux rgb
  622.  *  CLR_XYZ       rgb_aux[4]    (in)  - the color space definition,
  623.  *                                      3 primaries and white
  624.  *
  625.  * Creates the transformations between RGB color space the CLR_
  626.  *  routines are initialized for and another RGB color space as
  627.  *  specified in 'aux_rgb'.  The 'to' and 'from' matrices are
  628.  *  filled with the resulting transformations.  TRUE is returned
  629.  *  if successful, FALSE is returned if the CLR_ routines are not
  630.  *  initialized or if there is a singularity detected when the
  631.  *  transformation is being generated.
  632.  *
  633.  * The technique used is described in Sect 3.2, Colorimetry and
  634.  *  the RGB Monitor.
  635.  */
  636. int CLR_rgb_to_aux_rgb(to,from,rgb_aux)
  637. double  to[][3], from[][3];
  638. CLR_XYZ rgb_aux[];
  639. {   CLR_XYZ     rgb_tmp[4];
  640.     double      rgb_xyz_aux[3][3], xyz_rgb_aux[3][3];
  641.     int         color;
  642.     if (!init) return FALSE;
  643.     /* normalize the chromaticities of the auxiliary primaries
  644.      *  and white point.
  645.      */
  646.     for (color=RED; color<=WHITE; color++) {
  647.     rgb_tmp[color] = rgb_aux[color];
  648.     rgb_tmp[color].z = 1.0 - rgb_aux[color].x - rgb_aux[color].y;
  649.     }
  650.     /* get the transform between XYZ and the auxiliary RGB */
  651.     if (!clr__cspace_to_xyz(rgb_tmp, rgb_xyz_aux)) return FALSE;
  652.     if (!CLR_t_inverse(rgb_xyz_aux,xyz_rgb_aux)) return FALSE;
  653.     /* concatenate with the transforms for the RGB color space
  654.      *  that the CLR_ routine set is initialized to
  655.      */
  656.     (void)CLR_t_concat(xyz_rgb_aux, rgb_xyz_mat, to);
  657.     (void)CLR_t_concat(xyz_rgb_mat, rgb_xyz_aux, from);
  658.     return TRUE;
  659. }
  660.  
  661. /* ****************************************************************
  662.  * CLR_LAB CLR_xyz_to_lab (xyz)
  663.  *  CLR_XYZ     xyz     (in) - xyz color (.01<=Y<=1)
  664.  *
  665.  * Returns L*a*b* given XYZ.  Returns (0,0,0) if the CLR_ routines
  666.  *  are not initialized.  This transformation is described in
  667.  *  Section 3.3, Alternate Color Representations, Eq.(\*(sC),
  668.  *  and is taken from Judd and Wyszecki (1975) pg.320.  The
  669.  *  transformation is scaled for (.01 < Y < 1) for consistency
  670.  *  within the CLR_ routine set.
  671.  */
  672. CLR_LAB CLR_xyz_to_lab(xyz)
  673. CLR_XYZ     xyz;
  674. {   CLR_LAB     lab;
  675.     if (!init) {lab.l = lab.a = lab.b = 0.0; return lab;}
  676.     xyz.x *= 100.0;
  677.     xyz.y *= 100.0;
  678.     xyz.z *= 100.0;
  679.     lab.l = 25.0 * pow(((100.0*xyz.y)/white.y),.33333) - 16.0;
  680.     lab.a = 500.0 *
  681.     (pow(xyz.x/white.x,.33333) - pow(xyz.y/white.y,.33333));
  682.     lab.b = 200.0 *
  683.     (pow(xyz.y/white.y,.33333) - pow(xyz.z/white.z,.33333));
  684.     return lab;
  685. }
  686.  
  687. /* ****************************************************************
  688.  * CLR_LUV CLR_xyz_to_luv (xyz)
  689.  *  CLR_XYZ     xyz     (in) - xyz color (.01<=Y<=1)
  690.  *
  691.  * Returns L*u*v* given XYZ.  Returns (0,0,0) if the CLR_ routines
  692.  *  are not initialized.  This transformation is described in
  693.  *  Section 3.3, Alternate Color Representations, Eq.(\*(sD),
  694.  *  and is taken from Judd and Wyszecki (1975) pg.328.  The
  695.  *  transformation is scaled for (.01 < Y < 1) for consistency
  696.  *  within the CLR_ routine set.
  697.  */
  698. CLR_LUV CLR_xyz_to_luv(xyz)
  699. CLR_XYZ     xyz;
  700. {   CLR_LUV     luv;
  701.     double      u, v;
  702.     if (!init) {luv.l = luv.u = luv.v = 0.0; return luv;}
  703.     xyz.x *= 100.0;
  704.     xyz.y *= 100.0;
  705.     xyz.z *= 100.0;
  706.     luv.l = 25.0 * pow(((100.0*xyz.y)/white.y),.33333) - 16.0;
  707.     u = (4.0 * xyz.x) / (xyz.x + (15.0 * xyz.y) + (3.0 * xyz.z));
  708.     v = (9.0 * xyz.y) / (xyz.x + (15.0 * xyz.y) + (3.0 * xyz.z));
  709.     luv.u = 13.0 * luv.l * (u - ref_luv.u);
  710.     luv.v = 13.0 * luv.l * (v - ref_luv.v);
  711.     return luv;
  712. }
  713.  
  714. /* ****************************************************************
  715.  * CLR_t_concat (m1, m2, m3)
  716.  *  double        m1[3][3]      (in)  - matrix
  717.  *  double        m2[3][3]      (in)  - matrix to concat
  718.  *  double        m3[3][3]      (in)  - concatenated matrix
  719.  *
  720.  * Concatenate 'm1' to 'm2' resulting in 'm3'. In use, suppose you
  721.  *  have an XYZ to RGB and an RGB to YIQ matrix.  Concatenate the
  722.  *  RGB to YIQ matrix to the XYZ to RGB matrix to get an XYZ to
  723.  *  YIQ matrix.  Returns TRUE.
  724.  */
  725. int CLR_t_concat(m1,m2,m3)
  726. double  m1[][3], m2[][3], m3[][3];
  727. {   double  t1[3][3], t2[3][3];
  728.     LOAD_MAT(t1,m1);
  729.     LOAD_MAT(t2,m2);
  730.     {   int ii,jj;
  731.     for (ii=0; ii<=2; ii++)
  732.         for (jj=0; jj<=2; jj++)
  733.         m3[ii][jj] = (t1[ii][0] * t2[0][jj]) +
  734.                  (t1[ii][1] * t2[1][jj]) +
  735.                  (t1[ii][2] * t2[2][jj]);
  736.     }
  737.     return TRUE;
  738. }
  739.  
  740. /* ****************************************************************
  741.  * CLR_t_inverse (mat, inv_mat)
  742.  *  double        mat[3][3]     (in)  - matrix to be inverted
  743.  *  double        inv_mat[3][3] (mod) - inverted matrix
  744.  *
  745.  * Inverts 'mat' using Gaussian elimination.  Returns TRUE if
  746.  *  successful and FALSE if there is a singularity.
  747.  */
  748. int CLR_t_inverse (mat, inv_mat)
  749. double  mat[][3];
  750. double  inv_mat[][3];
  751. {   int     ii, jj, kk;
  752.     double  tmp_mat[3][3], tmp_d;
  753.  
  754.     for (ii=0; ii<=2; ii++)
  755.     for (jj=0; jj<=2; jj++) {
  756.         tmp_mat[ii][jj] = mat[ii][jj];
  757.         inv_mat[ii][jj] = (ii==jj ? 1.0 : 0.0);
  758.     }
  759.  
  760.     for (ii=0; ii<=2; ii++) {
  761.     for (jj=ii+1, kk=ii; jj<=2; jj++)
  762.         if (fabs(tmp_mat[jj][ii]) > fabs(tmp_mat[kk][ii]))
  763.         kk = jj;
  764.  
  765.     /* check for singularity */
  766.     if (tmp_mat[kk][ii] == 0.0) return FALSE;
  767.  
  768.     /* pivot - switch rows kk and ii */
  769.     if (kk != ii)
  770.         for (jj=0; jj<=2; jj++) {
  771.         tmp_d = tmp_mat[ii][jj];
  772.         tmp_mat[ii][jj] = tmp_mat[kk][jj];
  773.         tmp_mat[kk][jj] = tmp_d;
  774.         tmp_d = inv_mat[ii][jj];
  775.         inv_mat[ii][jj] = inv_mat[kk][jj];
  776.         inv_mat[kk][jj] = tmp_d;
  777.         }
  778.  
  779.     /* normalize the row - make the diagonal 1 */
  780.     for (tmp_d = 1.0 / tmp_mat[ii][ii], jj=0; jj<=2; jj++) {
  781.         tmp_mat[ii][jj] *= tmp_d;
  782.         inv_mat[ii][jj] *= tmp_d;
  783.     }
  784.  
  785.     /* zero the non-diagonal terms in this column */
  786.     for (jj=0; jj<=2; jj++)
  787.         if (jj != ii)
  788.         for (tmp_d = -tmp_mat[jj][ii], kk=0; kk<=2; kk++) {
  789.             tmp_mat[jj][kk] += tmp_mat[ii][kk] * tmp_d;
  790.             inv_mat[jj][kk] += inv_mat[ii][kk] * tmp_d;
  791.         }
  792.     }
  793.  
  794.     return TRUE;
  795. }
  796.  
  797. /* ****************************************************************
  798.  * clr__cspace_to_xyz (cspace, t_mat)
  799.  *  CLR_XYZ       cspace[4]   (in)  - the color space definition,
  800.  *                                      3 primaries and white
  801.  *  double        t_mat[3][3] (mod) - the color transformation
  802.  *
  803.  * Builds the transformation from a set of primaries to the CIEXYZ
  804.  *  color space.  This is the basis for the generation of the color
  805.  *  transformations in the CLR_ routine set.  The method used is
  806.  *  that detailed in Sect 3.2 Colorimetry and the RGB monitor.
  807.  *  Returns FALSE if there is a singularity.
  808.  */
  809. static int clr__cspace_to_xyz (cspace, t_mat)
  810. CLR_XYZ     cspace[];
  811. double      t_mat[][3];
  812. {   int     ii, jj, kk, tmp_i, ind[3];
  813.     double  mult, white[3], scale[3];
  814.  
  815.     /* normalize the white point to Y=1 */
  816.     if (cspace[WHITE].y <= 0.0) return FALSE;
  817.     white[0] = cspace[WHITE].x / cspace[WHITE].y;
  818.     white[1] = 1.0;
  819.     white[2] = cspace[WHITE].z / cspace[WHITE].y;
  820.  
  821.     for (ii=0; ii<=2; ii++) {
  822.     t_mat[0][ii] = cspace[ii].x;
  823.     t_mat[1][ii] = cspace[ii].y;
  824.     t_mat[2][ii] = cspace[ii].z;
  825.     ind[ii] = ii;
  826.     }
  827.  
  828.     /* gaussian elimination  with partial pivoting */
  829.     for (ii=0; ii<2; ii++) {
  830.     for (jj=ii+1; jj<=2; jj++)
  831.         if (fabs(t_mat[ind[jj]][ii]) >
  832.             fabs(t_mat[ind[ii]][ii])) {
  833.         tmp_i=ind[jj]; ind[jj]=ind[ii]; ind[ii]=tmp_i;
  834.         }
  835.     if (t_mat[ind[ii]][ii] == 0.0) return FALSE;
  836.  
  837.     for (jj=ii+1; jj<=2; jj++) {
  838.         mult = t_mat[ind[jj]][ii] / t_mat[ind[ii]][ii];
  839.         for (kk=ii+1; kk<=2; kk++)
  840.         t_mat[ind[jj]][kk] -= t_mat[ind[ii]][kk] * mult;
  841.         white[ind[jj]] -= white[ind[ii]] * mult;
  842.     }
  843.     }
  844.     if (t_mat[ind[2]][2] == 0.0) return FALSE;
  845.  
  846.     /* back substitution to solve for scale */
  847.     scale[ind[2]] = white[ind[2]] / t_mat[ind[2]][2];
  848.     scale[ind[1]] = (white[ind[1]] - (t_mat[ind[1]][2] *
  849.     scale[ind[2]])) / t_mat[ind[1]][1];
  850.     scale[ind[0]] = (white[ind[0]] - (t_mat[ind[0]][1] *
  851.     scale[ind[1]]) - (t_mat[ind[0]][2] * scale[ind[2]])) /
  852.     t_mat[ind[0]][0];
  853.  
  854.     /* build matrix */
  855.     for (ii=0; ii<=2; ii++) {
  856.     t_mat[0][ii] = cspace[ii].x * scale[ii];
  857.     t_mat[1][ii] = cspace[ii].y * scale[ii];
  858.     t_mat[2][ii] = cspace[ii].z * scale[ii];
  859.     }
  860.  
  861.     return TRUE;
  862. }
  863.  
  864. /* ****************************************************************
  865.  * CLR_exit()
  866.  *
  867.  * Completes use of the CLR_ routine set and frees any
  868.  *   allocated space
  869.  */
  870. CLR_exit()
  871. {
  872.     if (!init) return TRUE;
  873.     (void)CLR_exit_samples();
  874.     if (X_tristim != NULL) free((char *)X_tristim);
  875.     if (Y_tristim != NULL) free((char *)Y_tristim);
  876.     if (Z_tristim != NULL) free((char *)Z_tristim);
  877.     if (work_curve != NULL) free((char *)work_curve);
  878.     X_tristim = Y_tristim = Z_tristim = work_curve = NULL;
  879.     init = FALSE;
  880.     return TRUE;
  881. }
  882. /* ************************************************************* */
  883.