home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / cprog / illum.lha / clr_sample.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-05-11  |  11.1 KB  |  359 lines
  1. /* ****************************************************************
  2.  *                         clr_sample.c
  3.  * ****************************************************************
  4.  *  MODULE PURPOSE:
  5.  *      This module contains the routines for spectral sampling.
  6.  *      The operations include defining the sampling space,
  7.  *      sampling spectral curves, and transformations from spectral
  8.  *      sample space to RGB or XYZ.
  9.  *
  10.  *  MODULE CONTENTS:
  11.  *      CLR_init_samples    - initialize spectral sampling
  12.  *      CLR_num_samples     - returns the number of samples
  13.  *      CLR_spect_to_sample - sample a spectral curve
  14.  *      CLR_get_sample_rgb  - get the sample to RGB matrix
  15.  *      CLR_get_sample_xyz  - get the sample to XYZ matrix
  16.  *      CLR_reconstruct     - reconstruct a spectral curve
  17.  *      CLR_exit_samples    - finish with spectral sampling
  18.  *
  19.  *  NOTES:
  20.  *      > The CLR_ routines must be initialized (CLR_init()) before
  21.  *          using any of the sampling routines.
  22.  *
  23.  *      > When the CLR_SAMPLE_HALL method is selected, sampling
  24.  *          uses abutting, non-overlapping, box sample and
  25.  *          reconstruction functions are described in Section
  26.  *          3.5.2 Spectral Sampling Approaches, and in Hall (1983).
  27.  *          This provides continuous sampling between the low bound
  28.  *          of the first sample and the high bound of the last
  29.  *          sample.  For applications requiring discrete isolated
  30.  *          samples, user modification of these routines
  31.  *          is required.
  32.  *
  33.  *      > When the CLR_SAMPLE_MEYER method is used, 4 samples are
  34.  *          used as described in Meyer (1988).
  35.  */
  36. #include <stdio.h>
  37. #include "clr.h"
  38.  
  39. static int      sample_type = -1;
  40. static int      samples = 0;
  41. static int      *bounds = NULL;
  42. static double   XYZ_to_ACC[3][3] = {{-0.0177,   1.0090, 0.0073},
  43.                     {-1.5370,   1.0821, 0.3209},
  44.                     { 0.1946,  -0.2045, 0.5264}};
  45. static double   ACC_to_XYZ[3][3];
  46. static double   samp_to_ACC[3][4] =
  47.             {{0.00000,  0.18892,    0.67493,    0.19253},
  48.              {0.00000,  0.00000,    0.31824,   -0.46008},
  49.              {0.54640,  0.00000,    0.00000,    0.00000}};
  50.  
  51. extern char     *malloc();
  52.  
  53. /* ****************************************************************
  54.  * CLR_init_samples (method, num_samples, sample_bounds)
  55.  *  int     method          (in) - sampling method:
  56.  *                              CLR_SAMPLE_MEYER    Meyer (1988)
  57.  *                              CLR_SAMPLE_HALL     Hall (1983)
  58.  *  int     num_samples     (in) - number of sample functions
  59.  *  int     sample_bounds[] (in) - boundaries of the sampling
  60.  *                              functions.  There must be
  61.  *                              num_samples+1 bounds arranged in
  62.  *                              ascending order in this array.
  63.  *
  64.  * For the CLR_SAMPLE_HALL method the bound wavelength is included
  65.  *  in the sampling function.  For example, using 3 samples with
  66.  *  bounds at (411, 491, 571, 651), the actual samples are 411-490,
  67.  *  491-570, and 571-650.
  68.  *
  69.  * The CLR_SAMPLE_MEYER method uses a prescribed sampling with 4
  70.  *  samples.  The num_samples and sample_bounds arguments are
  71.  *  ignored.
  72.  *
  73.  * Returns TRUE if successful, FALSE if sample bounds are not valid
  74.  *  or sampling is previously initialized to some other value.
  75.  */
  76. CLR_init_samples (method, num_samples, sample_bounds)
  77. int     method;
  78. int     num_samples;
  79. int     *sample_bounds;
  80. {   int     ct;
  81.  
  82.     if (method == CLR_SAMPLE_MEYER) {
  83.     samples = 4;
  84.     sample_type = method;
  85.     CLR_t_inverse (XYZ_to_ACC, ACC_to_XYZ);
  86.     }
  87.     else if (method == CLR_SAMPLE_HALL) {
  88.     /* There are two ways to do this, one is to save a
  89.      *  complete sampling curve for each sample.  The
  90.      *  other is to do the computation on the fly.  Here
  91.      *  it's done on the fly.  Only the bounds are saved.
  92.      */
  93.     if (num_samples <= 0) return FALSE;
  94.     if ((bounds = (int *)malloc((unsigned)(sizeof(int) *
  95.         (num_samples + 1)))) == NULL) goto error;
  96.  
  97.     bounds[0] = sample_bounds[0];
  98.     for (ct=0; ct<num_samples; ct++) {
  99.         if (sample_bounds[ct+1] <= sample_bounds[ct])
  100.         goto error;
  101.         bounds[ct+1] = sample_bounds[ct+1];
  102.     }
  103.  
  104.     samples = num_samples;
  105.     sample_type = method;
  106.     }
  107.     else {
  108.     goto error;
  109.     }
  110.     return TRUE;
  111.  
  112. error:
  113.     samples = 0;
  114.     if (bounds != NULL) free((char *)bounds);
  115.     return FALSE;
  116. }
  117.  
  118. /* ****************************************************************
  119.  * CLR_num_samples()
  120.  *
  121.  * Returns the number of samples for which sampling is initialized.
  122.  *  Returns 0 if sampling is not initialized.
  123.  */
  124. CLR_num_samples()
  125. {   return samples;
  126. }
  127.  
  128.  
  129. /* ****************************************************************
  130.  * CLR_spect_to_sample (spectral, sample)
  131.  *  double  *spectral   (in)  - spectral curve to be sampled
  132.  *  double  *sample     (mod) - array to receive the sampled
  133.  *                          values.
  134.  * Samples 'spectral' and loads the sample values into 'sample'.
  135.  *  Returns TRUE if successful and FALSE if CLR_ or the sampling
  136.  *  has not been initialized
  137.  */
  138. CLR_spect_to_sample (spectral, sample)
  139. double      *spectral;
  140. double      *sample;
  141. {   int     cur_wl, ct, max_wl, cur_samp;
  142.     double  cur_sum;
  143.  
  144.     if (samples <= 0) return FALSE;
  145.     if ((cur_wl = CLR_get_min_wl()) < 0) return FALSE;
  146.     max_wl = CLR_get_max_wl();
  147.  
  148.     if (sample_type == CLR_SAMPLE_MEYER) {
  149.     /* sample at 456.4nm, 490.9nm, 557.7nm, and 631.4nm
  150.      */
  151.     if ((cur_wl > 456) || (max_wl < 457)) *sample++ = 0.0;
  152.     else *sample++ = spectral[456-cur_wl] + (0.4 *
  153.         (spectral[457-cur_wl] - spectral[456-cur_wl]));
  154.  
  155.     if ((cur_wl > 490) || (max_wl < 491)) *sample++ = 0.0;
  156.     else *sample++ = spectral[490-cur_wl] + (0.9 *
  157.         (spectral[491-cur_wl] - spectral[490-cur_wl]));
  158.  
  159.     if ((cur_wl > 557) || (max_wl < 558)) *sample++ = 0.0;
  160.     else *sample++ = spectral[557-cur_wl] + (0.7 *
  161.         (spectral[558-cur_wl] - spectral[557-cur_wl]));
  162.  
  163.     if ((cur_wl > 631) || (max_wl < 632)) *sample++ = 0.0;
  164.     else *sample++ = spectral[631-cur_wl] + (0.4 *
  165.         (spectral[632-cur_wl] - spectral[631-cur_wl]));
  166.     }
  167.     else {
  168.     for ( ;cur_wl<bounds[0]; cur_wl++, spectral++) ;
  169.  
  170.     for (cur_samp=1, cur_sum=0.0, ct=0; cur_wl<=max_wl;
  171.         cur_wl++) {
  172.         while (cur_wl>=bounds[cur_samp]) {
  173.         if (ct == 0) *sample++ = 0.0;
  174.         else *sample++ = cur_sum / ct;
  175.         if (++cur_samp > samples) return TRUE;
  176.         cur_sum = 0.0;
  177.         ct = 0;
  178.         }
  179.         cur_sum += *spectral++;
  180.         ct++;
  181.     }
  182.     *sample = cur_sum / ct;
  183.     }
  184.     return TRUE;
  185. }
  186.  
  187. /* ****************************************************************
  188.  * CLR_get_sample_rgb (matrix)
  189.  *  double  *matrix     (mod) - matrix to be filled.
  190.  *
  191.  * Returns the matrix for conversion from the sampled space to the
  192.  *  RGB the CLR_ routines have been initialized for.  The matrix
  193.  *  is a 3 x num_samples matrix.
  194.  */
  195. CLR_get_sample_rgb(matrix)
  196. double      *matrix;
  197. {
  198.     double      *xyz = NULL;
  199.     double      xyz_rgb[3][3];
  200.     int         ct;
  201.  
  202.     /* get the XYZ matrix, then transform it into RGB
  203.      */
  204.     if ((xyz = (double *)malloc((unsigned)(3 *
  205.     sizeof(double) * samples))) == NULL) goto error;
  206.     if (!CLR_get_sample_xyz(xyz)) goto error;
  207.     if (!CLR_get_xyz_rgb(xyz_rgb)) goto error;
  208.  
  209.     for (ct=0; ct<samples; ct++, matrix++, xyz++) {
  210.     matrix[0] = (xyz_rgb[0][0] * xyz[0]) +
  211.             (xyz_rgb[0][1] * xyz[samples]) +
  212.             (xyz_rgb[0][2] * xyz[2*samples]);
  213.     matrix[samples] = (xyz_rgb[1][0] * xyz[0]) +
  214.             (xyz_rgb[1][1] * xyz[samples]) +
  215.             (xyz_rgb[1][2] * xyz[2*samples]);
  216.     matrix[2*samples] = (xyz_rgb[2][0] * xyz[0]) +
  217.             (xyz_rgb[2][1] * xyz[samples]) +
  218.             (xyz_rgb[2][2] * xyz[2*samples]);
  219.     }
  220.  
  221.     free((char *)xyz);
  222.     return TRUE;
  223.  
  224. error:
  225.     if (xyz != NULL) free((char *)xyz);
  226.     return FALSE;
  227. }
  228.  
  229. /* ****************************************************************
  230.  * CLR_get_sample_xyz (matrix)
  231.  *  double  *matrix     (mod) - matrix to be filled.
  232.  *
  233.  * Returns the matrix for conversion from the sampled space to
  234.  *  CIEXYZ. The matrix is a 3 x num_samples matrix.
  235.  */
  236. CLR_get_sample_xyz(matrix)
  237. double      *matrix;
  238. {
  239.     int     min_wl, max_wl, cur_wl, cur_samp;
  240.     double  *reconst, *cur_r, fill, *samp_mat;
  241.     CLR_XYZ xyz;
  242.  
  243.     if (samples <= 0) return FALSE;
  244.     if (sample_type == CLR_SAMPLE_MEYER) {
  245.     /* concatenate the sample to ACC matrix with the ACC_to_XYZ
  246.      *  matrix.  The divide by 1.057863 is a normalization so
  247.      *  than an identity curve has a Y value of 1.0 following
  248.      *  the conventions used in the CLR_routines.
  249.      */
  250.     samp_mat = samp_to_ACC[0];
  251.     for (cur_samp=0; cur_samp<samples; cur_samp++,
  252.         matrix++, samp_mat++) {
  253.         matrix[0] = ((ACC_to_XYZ[0][0] * samp_mat[0]) +
  254.         (ACC_to_XYZ[0][1] * samp_mat[samples]) +
  255.         (ACC_to_XYZ[0][2] * samp_mat[2*samples])) /
  256.         1.057863;
  257.         matrix[samples] = ((ACC_to_XYZ[1][0] * samp_mat[0]) +
  258.         (ACC_to_XYZ[1][1] * samp_mat[samples]) +
  259.         (ACC_to_XYZ[1][2] * samp_mat[2*samples])) /
  260.         1.057863;
  261.         matrix[2*samples] = ((ACC_to_XYZ[2][0] * samp_mat[0]) +
  262.         (ACC_to_XYZ[2][1] * samp_mat[samples]) +
  263.         (ACC_to_XYZ[2][2] * samp_mat[2*samples])) /
  264.         1.057863;
  265.     }
  266.     }
  267.     else {
  268.     min_wl = CLR_get_min_wl();
  269.     max_wl = CLR_get_max_wl();
  270.     /* allocate space for the reconstruction function
  271.      */
  272.     if ((reconst = (double *)malloc(
  273.         (unsigned)(sizeof(double) *
  274.         (max_wl - min_wl + 1)))) == NULL) goto error;
  275.  
  276.     /* Build each reconstruction function and sample it into xyz.
  277.      *  Load the values into the matrix
  278.      */
  279.     for (cur_samp=0; cur_samp<samples; cur_samp++, matrix++) {
  280.         cur_r = reconst;
  281.  
  282.         if (cur_samp == 0) fill = 1.0;
  283.         else fill = 0.0;
  284.  
  285.         for (cur_wl=min_wl ;(cur_wl<bounds[cur_samp]) &&
  286.         (cur_wl<max_wl); cur_wl++) *cur_r++ = fill;
  287.         for ( ;(cur_wl<bounds[cur_samp+1]) &&
  288.         (cur_wl<=max_wl); cur_wl++) *cur_r++ = 1.0;
  289.         if (cur_samp == (samples-1)) fill = 1.0;
  290.         else fill = 0.0;
  291.         for ( ;cur_wl<=max_wl; cur_wl++) *cur_r++ = fill;
  292.         xyz = CLR_spect_to_xyz (reconst);
  293.         matrix[0] = xyz.x;
  294.         matrix[samples] = xyz.y;
  295.         matrix[2*samples] = xyz.z;
  296.     }
  297.     }
  298.  
  299.     free((char *)reconst);
  300.     return TRUE;
  301.  
  302. error:
  303.     if (reconst != NULL) free((char *)reconst);
  304.     return FALSE;
  305. }
  306.  
  307. /* ****************************************************************
  308.  * CLR_reconstruct (sample, spectral)
  309.  *  double  *sample     (in)  - the sample values
  310.  *  double  *spectral   (mod) - the reconstructed spectral curve
  311.  *
  312.  * Reconstructs a spectral curve from the sample values.  The
  313.  *  reconstruction functions are box functions resulting a a step
  314.  *  function as the reconstructed curve.
  315.  */
  316. CLR_reconstruct (sample, spectral)
  317. double      *sample, *spectral;
  318. {   int     cur_wl, max_wl, cur_sample;
  319.     if (samples <= 0) return FALSE;
  320.     if (sample_type == CLR_SAMPLE_MEYER) {
  321.     return FALSE;
  322.     }
  323.     else if (sample_type == CLR_SAMPLE_HALL) {
  324.     cur_wl = CLR_get_min_wl();
  325.     max_wl = CLR_get_max_wl();
  326.     cur_sample = 1;
  327.  
  328.     while (cur_wl<=max_wl) {
  329.         if (cur_wl < bounds[cur_sample])
  330.         *spectral++ = *sample;
  331.         else if (cur_sample >= samples)
  332.         *spectral++ = *sample;
  333.         else {
  334.         cur_sample++;
  335.         sample++;
  336.         *spectral++ = *sample;
  337.         }
  338.         cur_wl++;
  339.     }
  340.     }
  341.     return TRUE;
  342. }
  343.  
  344. /* ****************************************************************
  345.  * CLR_exit_samples()
  346.  *
  347.  * Complete use of spectral sampling, free any allocated space.
  348.  */
  349. CLR_exit_samples ()
  350. {   sample_type = -1;
  351.     samples = 0;
  352.     if (bounds != NULL) {
  353.     free((char *)bounds);
  354.     bounds = NULL;
  355.     }
  356.     return TRUE;
  357. }
  358. /* ************************************************************* */
  359.