Help for TRUCOLOR

PURPOSE
"trucolor" transforms
images taken through several filters into a color image in which an
attempt is made to reproduce accurately certain designated spectra ("special
colors").  At this time, the only output devices supported by this program are
the CONRAC 7211 color monitor, and the MDA Color Fire film recorder. 


EXECUTION
  The program is executed by specifying up to 10 images
taken through different filters, up to 10 special
colors, and, three output files (red, green, and blue), 
a conversion factor for each input image, and an image display device.

  The input images should be 
reflectance images if the specified DEVICE is a film recorder, and radiance 
images if the specified DEVICE is a color monitor. 

  When the program begins it will look in your local directory for
a file called "trucolor".TV or "trucolor".FILM . If it finds this file
it will read it to acquire the Yxy to DN interpolation table. If it
cannot find it it will recompute the table (taking several cpu minutes)
and then write it out for another execution of "trucolor" to find.
The file is about the size of a flight image.

  Also when the program begins it will look in your local directory
for a file called CUBE.TV or CUBE.FILM. If it finds such a file it will
read it and use the information to replace an internal table which
tabulates the Yxy values for 125 combinations of DN values, thus
modelling the device in question (either the tv or film recorder).
All values are real format in ascii.
If you create such a file it is ordered as 125 records with 3 real
values per record in the order Y (tristimulus value) x (chromaticity)
and y (chromaticity). 
The order of the sets of 3 points is:
dnred=0,64,128,192,255  for green=0  blue=0
         "                        64      0  
         "                        128     0  
         "                        192     0  
         "                        255     0  
then repeat the above for blue=64
then repeat the above for blue=128
then repeat the above for blue=192
then repeat the above for blue=255

  If you provide four output files the last output will contain
a chromaticity map histogram of all the data points in the inputs.

OUTPUT FILES:
All outputs are in BYTE format.

  If, for any reason, the output dn value cannot be computed it will be
set to red=green=blue=0(zero). This is true for ANY saturation in ANY
part of a red,green,blue value pixel, for Yxy values outside the
device range, or negative Y.

  The fourth output (optional) contains a chromaticity 2-d histogram.
The file is 256 by 256 byte format on a 100 dn background. 
The x chromaticity axis is in the sample direction.
The y chromaticity axis is in the line direction.
Each coordinate value for the display device mapping dn values into
Yxy values is recorded as ZERO dn or black. There are 125 points which
are black corresponding to the 125 values in the lookup table.
Each input dn value in the input pictures
is converted to Yxy and the (x,y) portion is accumulated above the
100 dn gray level as a 2-d histogram. This histogram will reveal the
clustering of colors in the images. Values outside the triangle formed
by the 125 black points will be set to zero dn in the output images
since they lie outside the available color range of the recording
device.

EXAMPLE
In the following examples the user creates output files, OUT.*, using input
images INP.*, and special colors, 11/7, 7/6, and 5/7.
The resulting images will be displayed on a
Conrac 7211 monitor. In this mode (radiance) Galsos must produce images
in radiance units.

  VICAR>trucolor INP=(INP.ORA,INP.CLR,INP.GRE) OUT=(OUT.RED,OUT.GRE,OUT.BLU)    VICAR>+ COLORS=(5,4,3) 'TV
 
WRITTEN BY: 		J Lorre 1/20/92

Made portable for UNIX  A Scop (CRI) 1/2/95

INTERNAL TABLES:
  You will note in the program two tables labelled:
Yxy_table  (for the film recorder)
Yxy_tv     (for the TV monitors)
These tables contain the Yxy values at each 64 dn levels for all
combinations of input dn values for blue,green, and red. You can
recreate these tables by measuring the Yxy values from the
photographic paper or from the TV monitor. The following .pdf
files were created to assist in this process.

TV monitor procedure:

procedure
refgbl $echo
body
let _onfail="continue"
let $echo="yes"
local (red,green,blue) integer
local (epa,return) string
write "specify display device: example EPA2"
getpar epa 
use &epa
vids
jcolor
jerase
for red=0,64,128,192,255
  gen &"red".img nl=200 ns=200 ival=&red linc=0 sinc=0
end-for
write "Center the probe on this patch"
jdisp (255.img,255.img,255.img) loc=(200,200)
getpar return
write "Now we cycle through all 125 combinations"
for blue=0,64,128,192,255
  for green=0,64,128,192,255
    for red=0,64,128,192,255
      write "red_dn=&red green_dn=&green blue_dn=&blue"
      jdisp (&"red".img,&"green".img,&"blue".img) loc=(200,200)
      getpar return
    end-for
  end-for
end-for
for red=0,64,128,192,255
  dcl del &"red".img;1
end-for
jstop
free &epa
end-proc

The above procedure displays patches at the center of a 512 TV monitor
in the order input to the tables. You may tape the Minolta probe to
the TV tube and write the Yxy values down as they come up.

Film procedure:

procedure
refgbl $echo
body
let _onfail="continue"
let $echo="yes"
!
! creates the FILM image to calibrate the film recorder
gen out=0.img nl=500 ns=100 ival=0 linc=0 sinc=0
gen out=1.img nl=500 ns=100 ival=64 linc=0 sinc=0
gen out=2.img nl=500 ns=100 ival=128 linc=0 sinc=0
gen out=3.img nl=500 ns=100 ival=192 linc=0 sinc=0
gen out=4.img nl=500 ns=100 ival=255 linc=0 sinc=0
fastmos inp=(0.img,1.img,2.img,3.img,4.img) out=r.img    nl=500 ns=500    off1=(1,1) off2=(1,101) off3=(1,201) off4=(1,301) off5=(1,401)
fastmos inp=(r.img,r.img,r.img,r.img,r.img) out=v2$scratch:red.img    nl=1500 ns=1000    off1=(1,1) off2=(1,501) off3=(501,1) off4=(501,501) off5=(1001,1)
flot inp=r.img out=g.img 'clock
fastmos inp=(g.img,g.img,g.img,g.img,g.img) out=v2$scratch:green.img    nl=1500 ns=1000    off1=(1,1) off2=(1,501) off3=(501,1) off4=(501,501) off5=(1001,1)
gen out=0.img nl=500 ns=500 ival=0 linc=0 sinc=0
gen out=1.img nl=500 ns=500 ival=64 linc=0 sinc=0
gen out=2.img nl=500 ns=500 ival=128 linc=0 sinc=0
gen out=3.img nl=500 ns=500 ival=192 linc=0 sinc=0
gen out=4.img nl=500 ns=500 ival=255 linc=0 sinc=0
fastmos inp=(0.img,1.img,2.img,3.img,4.img) out=v2$scratch:blue.img    nl=1500 ns=1000    off1=(1,1) off2=(1,501) off3=(501,1) off4=(501,501) off5=(1001,1)
barne inp=("v2$scratch:red.img,v2$scratch:green.img,    v2$scratch:blue.img") miplbox=21 primary=42057
dcl del %.img;1
!
end-proc

The above procedure creates a single print for the Minolta chromaticity
meter held by Dave Deats at the photo lab. The print is ordered as 5
squares like this:

*******************
*        *        *
*   0    *   64   *
*        *        *
*******************
*        *        *
*  128   *  192   *
*        *        *
*******************
*        *
*  255   *
*        *
**********

where the numbers are the BLUE dn values.
Each square has 25 patches arranged as a 5 by 5 grid with red increasing
from left to right 0 64 128 192 255, and green increasing from top down
0 64 128 192 255.


THEORY:
                          Color Calibration

   This document describes the process of color calibration for the
purpose of producing color products which look the same to the eye
as the original. 

   We first address the camera signal (dn) and the computation of U 
and I/F which are the decalibrated camera response. The response (dn) 
of a ccd detector is:

    t*O*a
dn= ----- integral(S*T*F*U*P*dl) + dn0                        (1)
     g*c
 
where:
dn=camera output number
g =gain state ratio factor
t =exposure time (sec)
c =high gain conversion factor (e/pixel/dn)
O =optics solid angle (sr)
a =area of ccd pixel (cm**2/pixel)
S =ccd spectral sensitivity (e/photon)
T =optics spectral transmission
F =filter spectral transmission
U =spectral scene radiance (watts/cm**2/sr/nm)
P =wavelength (nm) /1.9862e-16 (photons/watt-sec)
l =wavelength (nm)
dn0=dark current output number

   We wish to decalibrate this sensor before using the response (dn).
To do this we must extract the radiance U from the integral and solve
for it. The only way to extract U is to assume it is independent of 
wavelength. This is a serious error but one which must be made if
only one filter position is available. Typically we have three or more
filters if we are to perform color reconstruction, but for now we will
ignore this problem.

   Solving (1) for U we have:

      (dn-dn0)*g*c
U=-------------------------         watts/cm**2/sr/nm         (2)
  t*O*a*integral(S*T*F*P*dl)

or, in terms of reflectance I/F we have:

      (dn-dn0)*g*c*pi
I/F=-----------------------------     dimensionless albedo     (3)
     t*O*a*integral(S*T*F*H*P*dl)

where H is the solar irradiance incident on the surface.

To simplify things we denote from now on in the text both U
and I/F as symbol Q. In general U is used to produce color
for an additive device like a color TV monitor and I/F is used
for a subtractive device like a film recorder. 

We adopt a set of i special colors with known spectra C(i).
These spectra represent colors which we intend the camera to
reproduce with some precision. For example Kodak uses sky blue,
white cloud, green grass, and skin color as it's four special
colors for film process control. We could use Jupiter 
characteristic spectra.
The tristimulus values of these colors are computed from:

X(i)=integral(C(i)*x*dl)
Y(i)=integral(C(i)*y*dl)     1<=i<=I                           (4)
Z(i)=integral(C(i)*z*dl)

where: x,y,z are the standard CIE color matching functions.

We expose the camera to each of the colors to record the Q(ij) 
responses given by equation(2 or 3) in each of J filter positions.

Because the form of equations (1) and (4) are similar and linear it
is possible to represent the tristimulus values as a linear
combination of the camera responses Q(ij) to each special color i
in filter position j:

X(i)=A1*Q(i1)+A2*Q(i2)+A3*Q(i3)...AJ*Q(iJ)    
Y(i)=B1*Q(i1)+B2*Q(i2)+B3*Q(i3)...BJ*Q(iJ)   1<=i<=I        (5)
Z(i)=C1*Q(i1)+C2*Q(i2)+C3*Q(i3)...CJ*Q(iJ)    

There are as many X equations as there are special colors and as many
A terms as filters. We can solve for the A,B,C coefficients from each
group of the above equations. For input Q(ij) values taken through the
same filters used to solve equations (5) we can compute the tristimulus
values X,Y,Z. These are the values the eye should see.

To calibrate the display device we create a numerical step target
with discrete combinations of red,green,blue spanning the range of
numeric values. The target is displayed on the TV monitors and on
a print. For each discreet step we compute the tristimulus values:

for TV
X=integral(U*x*dl)
Y=integral(U*y*dl)                                   (6)
Z=integral(U*z*dl)

for prints
X=integral(R*P*x*dl)
Y=integral(R*P*y*dl)                                 (7)
Z=integral(R*P*z*dl)

where:
U =the spectral radiance of the TV monitor.
R =the spectral reflectance of the print.
P =a standard illuminant source like D65.

We know by now the dn(blue,green,red) versus the X,Y,Z tabular
relationship for the output device. To facilitate interpolation it is
necessary to create a regularly spaced interpolation table which
returns output blue,green,red dn's from a desired input X,Y,Z.

To create a color product from a set of J input filters we do:

First:
1. Select a set of I representative special colors for which
   camera responses are known (from equations 2 or 3) and
   tristimulus values are known (from equation 4).
2. Solve equation (5) for the 3*J coefficients A(j),B(j),C(j).

Then, for each pixel:
3. From equation (5) compute the X,Y,Z tristimulus values.
4. Interpolate in the 3-d table of X,Y,Z values to obtain the
   output red,green,blue dn values.

   Scaling:

   Scaling of the results may be required in order to produce
unsaturated images. It will be evident from the interpolation table 
which gives the output dn versus tristimulus values what the
tristimulus value limit is (corresponding to dn's of 255). 
Simply scaling the tristimulus values such that the luminance
Y does not exceed a limit will not be satisfactory. This is 
because the maximum Y is a function of the chromaticity coordinates
(x,y). To perform scaling we must:

1. Convert the XYZ values to Yxy using the relation
      x=X/(X+Y+Z) and y=Y/(X+Y+Z)
   We perform this BOTH for input data and for the interpolation
   table.
   We retain Y as the luminance since it is defined as such
   if we use the CIE 1931 color matching functions.

2. For each point xy in the data search the interpolation
   table for the largest possible Y value. Call it Ymax.
   Compute the smallest ratio of Ymax/Y. Call it Sc.

3. When processing images convert XYZ to Yxy and impose
   scaling by multiplying each Y by Sc.

Note that I/F reflectance may not require scaling.

   Discussion

   We must present an argument for the validity of using U or I/F
since we have approximated these functions as independent of
wavelength. Equation (5) is justified because it is a linear
representation of responses which are themselves linear integrals.
Any approximation which is linear such as substituting a flat
function for the scene radiance can be compensated for by the
three sets of least squares fit coefficients A,B, and C.
The important criterion is that I/F be a measurable quantity.
Since dn is proportional with exposure and I/F is proportional
with dn then I/F is measurable and correlatable to X,Y, and Z.
Precisely the same method is used in determining chemical
abundances from spectra such as in NIMS. Equation (5) locks
the camera response to XYZ for the colors of interest.
   One way to get into trouble with equation (5) is to include
filter positions which fall partially out of the spectral
band of the color matching functions x,y,z. Any IR filter would
cause this problem. Another way is to include filters which 
overlap so much spectrally that the images they produce are
nearly identical. Both these cases can be mitigated to some
extent by using the first three terms in the principal 
component transformation of the input images. We must use the
same rotation matrix obtained from the special colors to rotate
the images of unknown scenes. 
   Equation (5) provides a least squares fit between 
the U or I/F values and tristimulus values, both obtained from
special color targets. Later on these same A,B,C coefficients
will be used to convert from U or I/F to tristimulus values for
each pixel in the image. It is imperative that identical
means be used to compute U or I/F (Q) in both cases. Thus the
cameras must be exposed to each of the special colors rather than
U or I/F being computed theoretically. The point is that both U
and I/F are wrong ! They are wrong because of the limitations imposed
by broad band radiometry. To get consistent results we must make the
same error at both the calibration and the evaluation ends.
   One can compute the I/F values theoretically for the special
colors instead of exposing the camera to them. As discussed this
is not correct but, for the sake of completeness, we present
the formula:

      integral(S*T*F*U*dl)
I/F =-------------------------
      integral(S*T*F*H*dl)

where H= The incident spectral irradiance to a surface.
      U= The reflected spectral radiance from the surface.
      U= H*R/pi, where R= The spectral reflectance of the surface.

   The color cube approach to calibrating the film recorder or the
TV monitor rather than modelling these devices as logarithmic or
exponential is, I think, necessary. An interpolation table can 
model any device AND can tell if the requested tristimulus values are
beyond those accessible for the device. Film recorders+film+paper
are highly nonlinear creations with color crosstalk between channels.

   Alternative color method.

   There is an alternative method to color reconstruction.
This is to estimate the true spectral distribution of U or
I/F from a synthesis of the U or I/F values obtained through a
set of different camera filters (the very case we have). In this
scenario we attempt to guess the spectral radiance U which 
solves equation (1) for the J dn values AT EACH PIXEL.
Then we can integrate this function using equation (6) to
achieve the tristimulus values at each pixel directly
without using equation (5). This has been done by
Steve Wall JGR,Vol 8,No 28,pp 4401-4411. In this case special
colors are used more as a means of refining the method of
estimating U (since it is known) and of verifying results.

   Required calibration data:

1. Special colors targets imaged by the camera.
   This includes knowledge of the illuminant and target reflectance
   or transmittance curves, and the camera I/F or U responses.

2. A matrix spanning all dn ranges for calibrating the film recorder
   and TV monitor. This includes computations of the tristimulus
   values at each matrix point. A 5 by 5 by 5 cube would be adequate.

PROGRAM HISTORY

Written by Jean Lorre circa 1990
Ported by CRI, Jan 2 1995
Current cognizant programer: Jean Lorre
Revision history:
  Sep 19 95  GMY  Replace call to xveaction with 'SA' parameters in xvopen
                  to avoid abending in xlget when 'EDRTAPE' item is missing
                  from VICAR label (FR 89820)