Meet iCAM: An Image Color Appearance Model
Garrett M. Johnson
Mark D. Fairchild

Tone Reproduction for High Dynamic Range Images

Input

First read in the input image, and display the RGB version. We first need to do a little Mathematica house cleaning to get the path of the input image.

nbinfo = NotebookInformation[EvaluationNotebook[]]; 
    dir = ("FileName" /. nbinfo /. 
          FrontEnd`FileName[d_List, nam_, ___] :> 
            ToFileName[d]);

Read in the High-dynamic range image. This image was originally generated by Debevec, and can be found at 		http://www.debevec.org/Research/HDR/

The Images are originally in Greg Larsons' Radiance image format. Radiance was used to convert these images to a standard binary image format, and then embedded in an HDF format. Radiance can be found at:
 		http://radsite.lbl.gov/radiance/HOME.html
 information on the HDF image format can be found at: 
 		http://hdf.ncsa.uiuc.edu/
 
 The HDF image is in reverse order from what Mathematica expects, so we must reverse the input of the Import function.

rgbImage = 
        Reverse[Import["memorial_small.hdf",  
            Path -> dir]];

Render the HDR linear image. Except for a few highlights, it is almost completely dark.

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
(rgbImage/Max[rgbImage]), {2}]]],  
        AspectRatio -> Automatic];

Find the dimensions of the image, for future calculations

imDim = Dimensions[rgbImage];

imWidth = imDim[[1]];

imHeight  =  imDim[[2]];

Gamma Correction

Generally this stage involves linearizing the image, before transforming to XYZ values. Since the HDR image is already linear, this is not necessary. We will normalize the data to be between 0-1, and display it with a gamma correction, just to illustrate that  a gamma correction is not enough for HDR images.

γ = 1.8;

dcRGB = ((rgbImage/Max)[rgbImage]);

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
(dcRGB^((1/γ))), {2}]]],  AspectRatio -> Automatic];

RGB to XYZ Color Conversion

Convert the linearized RGB signal into XYZ using a standard sRGB matrix normalized for D65. Note that the Mathematica matrix convention is transposed from typical matrix descriptions.

sRGBmat=(({{0.4124, 0.3576, 0.1805}, {0.2126, 0.7152, 0.0722}, {0.0193, 0.1192, 0.9505}}));

xyzIm =  (dcRGB . Transpose[sRGBmat]);

Plot the XYZ images side-by-side for examination

Show[
        GraphicsArray[
{ListDensityPlot[
              xyzIm[[All, All, 1]], AspectRatio -> Automatic, 
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity], 
ListDensityPlot[
              xyzIm[[All, All, 
                  2]], 
AspectRatio -> Automatic,  
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity], 
ListDensityPlot[
              xyzIm[[All, All, 
                  3]], 
AspectRatio -> Automatic, 
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity]}], 
AspectRatio -> 
          Automatic, 
DisplayFunction -> $DisplayFunction];

test = 
        xyzIm[[All, All, 1]]/Max[xyzIm[[All, All, 1]]];

Isolate the individual X, Y, and Z channels

{xIm = Take[xyzIm[[All, All, 1]]], 
yIm = 
          Take[xyzIm[[All, All, 2]]], 
zIm = 
          Take[xyzIm[[All, All, 3]]]};

Max[yIm]

Min[yIm]

Determine Image "White" for Chromatic Adaptation

Create a smoothing kernel to use as the low-pass filter. First step is to load the image processing library, to allow for the discrete convolution. This can be replaced with frequency filtering, if desired.

<< ImageProcessing`

kernSize =  (imWidth/1);

Define a Gaussian convolution kernel

kern = 
        N[Table[(E^(-(((Sqrt[(x^2) + (y^2)])/((((kernSize/40)))^2))))), {x, -((kernSize - 
                      1)/2), ((kernSize - 1)/2)}, {y, -((kernSize - 
                      1)/2), ((kernSize - 1)/2)}]];

Normalize the convolution kernel so that it sums to 1.

kern = kern/Apply[Plus, Flatten[kern]];

yWhite = 
        DiscreteConvolve[yIm,  kern,  Centered -> True];

xWhite = yWhite;

zWhite = yWhite;

Plot the blurred XYZ images

Show[
        GraphicsArray[
{ListDensityPlot[xWhite, 
              AspectRatio -> Automatic, 
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity], 
ListDensityPlot[
              yWhite, 
AspectRatio -> Automatic,  
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity], 
ListDensityPlot[
              zWhite, 
AspectRatio -> Automatic, 
              Mesh -> False, 
Frame -> False,  
              DisplayFunction -> 
                Identity]}], 
AspectRatio -> 
          Automatic, 
DisplayFunction -> $DisplayFunction];

Combine the individual XYZ images back into a single [x,y,3] image.

xyzWhite = xyzIm;
xyzWhite[[All, All, 1]] = xWhite ;  
    xyzWhite[[All, All, 2]] = yWhite ; 
    xyzWhite[[All, All, 3]] = zWhite ;

Chromatic Adaptation

Set up some variables for the chromatic adaptation transform. This utilizes a slightly modified version of the CIECAM02 transform to include non-linear exponential functions, if desired.

Adapting Luminance

La = 100.0 ;

Surround Conditions

F = 1.0 ;

Degree of Adaptation

Dfact = 
      F (1 - (1/(1 + 2 (La^((1.0/4.0))) + ((La^2)/300.0))))

CIEDf = 
      F (1 - (1/
                  3.6) (E^((((-La - 42)/92)))))

Df =  0.1

Adapting Whitepoint

whitePoint=(({{95.05, 100.0, 108.9}}))/100;

Chromatic Adaptation Matrix

mCAT02=(({{0.7328, 0.4296, -0.1624}, {-0.7036, 1.6974, 0.0061}, {0.0030, 0.0136, 0.9834}}));

iCAT02 = Inverse[mCAT02];

MatrixForm[iCAT02]

Calculate RGB values for input image, blurred image, and whitepoint

rgbI = xyzIm . Transpose[mCAT02];

rgbWhitePoint = whitePoint . Transpose[mCAT02]

Calculate the RGB of the blurred "whitepoint" image. First define a function to threshold the minimum of the image to 1. This means the darkest RGB adapting value is 1, and prevents division by zero, or division by a small number.

threshold[x_] := If[x >  .001, x,  .001];

rgbWhite = xyzWhite . Transpose[mCAT02];

rgbWhite = Max[rgbI] rgbWhite/Max[rgbWhite] ;

Calculate the Maximum R, G, B value

{maxR = Max[rgbI[[All, All, 1]]],  
      maxG = Max[rgbI[[All, All, 2]]], 
maxB = 
        Max[rgbI[[All, All, 3]]]}

Calculate Adapted RGB Signals

Rc = (((Df rgbWhitePoint[[1, 1]])/rgbWhite)[[All, All, 1]] + 1 - Df) *rgbI[[All, All, 1]];

Gc = (((Df rgbWhitePoint[[1, 2]])/rgbWhite)[[All, All, 2]] + 1 - Df) *rgbI[[All, All, 2]];

Bc = (((Df rgbWhitePoint[[1, 3]])/rgbWhite)[[All, All, 3]] + 1 - Df) *rgbI[[All, All, 3]];

Convert Adapted RGB Signals back to XYZ, Examine the RGB color image for reference

RGBc = rgbWhite * 0 ;

{RGBc[[All, All, 1]] = Rc,  RGBc[[All, All, 2]] = Gc,  
        RGBc[[All, All, 3]] = Bc};

xyzAdapt = RGBc . Transpose[iCAT02];

rgbAdapt = xyzAdapt . Transpose[Inverse[sRGBmat]];

rgbAdapt = Map[threshold,  rgbAdapt,  {3}];

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
(((rgbAdapt/Max[rgbAdapt]))^((1/1.8))), {2}]]],  
        AspectRatio -> Automatic];

Convert to IPT color space

The IPT matrix is defined to transform LMS cone signals into the opponent IPT space. The LMS functions are modified using a non-linear exponent function. Typically this is a value of 0.43, but in iCAM we modify this exponent based on a low-passed Y image.

lmsIPTm=(({{0.4000, 0.4000, 0.2000}, {4.4550, -4.8510, 0.3960}, {0.8056, 0.3572, -1.1628}}));

xyzLMSm=(({{0.3971, 0.68898, -0.07868}, {-0.22981, 1.18340, 0.04641}, {0.0, 0.0, 1.0}}));

Now convert the XYZ adapted image into LMS cone space

lmsImage = xyzAdapt . Transpose[xyzLMSm];

Determine the exponents , based on the luminance of the input image

iptKernSize = (imWidth/1);

iptKern = 
        N[Table[(E^(-((Sqrt[(x^2) + (y^2)])/((((iptKernSize/40)))^2)))), {x, -((iptKernSize - 
                      1)/2), ((iptKernSize - 1)/2)}, {y, -((iptKernSize - 1)/2), ((iptKernSize - 1)/2)}]];

iptKern = iptKern/Apply[Plus,  Flatten[iptKern]];

yLow = 
        DiscreteConvolve[yIm,  iptKern,  
          Centered -> True];

ListDensityPlot[yLow,  Frame -> False,  Mesh -> False,  
        AspectRatio -> Automatic];

Efact[
        L_] = ((1.0/1.7)) (0.2 ((((1/(5 L  + 1.0))))^4) (5 L) +  0.1 (((1 - ((((1/(5 L  + 1.0))))^4)))^2) (((5 L))^(1/3)))

iptExp = Efact[(1000 yLow)];

Create a minimum threshold for the exponential, otherwise we will clip dark values to zero.

Max[iptExp]* .43

minThresh[x_] := If[x < 0.3,  0.3,  x];

iptExp = Map[minThresh, iptExp,  {2}];

lmsImNL = (lmsImage^(iptExp* .43));

Convert the non-linear LMS signals into IPT space

iptImage = lmsImNL . Transpose[lmsIPTm];

Plot the "I" dimension, along with the Luminance (Y) image

Show[
        GraphicsArray[{ListDensityPlot[
              iptImage[[All, All, 1]]^((1/1)),  Frame -> False,  
              Mesh -> False,  AspectRatio -> Automatic,  
              DisplayFunction -> Identity],  
            ListDensityPlot[xyzIm[[All, All, 2]],  
              Frame -> False,  Mesh -> False,  
              AspectRatio -> Automatic,  
              DisplayFunction -> Identity]}], 
        DisplayFunction -> $DisplayFunction];

Inverse IPT transform

iptLMSm = Inverse[lmsIPTm];

lmsXYZm = Inverse[xyzLMSm];

xyzRGBm = Inverse[sRGBmat];

MatrixForm[iptLMSm]

MatrixForm[lmsXYZm]

MatrixForm[xyzRGBm]

lmsIPTimage = iptImage . Transpose[iptLMSm];

lmsIPTimage = Map[threshold,  lmsIPTimage,  {3}];

xyzIPTimage = 
        (lmsIPTimage^((1/0.43))) . Transpose[lmsXYZm];

Inverse Chromatic Adaptation

rgbIPTim = xyzIPTimage . Transpose[mCAT02];

monWhite=(({{1.0, 1.0, 1.0}})) ;

rgbWhite = monWhite . Transpose[mCAT02]

Note we are setting the Dfactor to one for the inverse chromatic adaptation transform, so that we are completely adapting.

Df = 1.0

Rc = (((Df rgbWhite[[1, 1]])/rgbWhitePoint)[[1, 1]] + 1 - 
              Df) *rgbIPTim[[All, All, 1]];

Gc = (((Df rgbWhite[[1, 2]])/rgbWhitePoint)[[1, 2]] + 1 - 
              Df) *rgbIPTim[[All, All, 2]];

Bc = (((Df rgbWhite[[1, 3]])/rgbWhitePoint)[[1, 3]] + 1 - 
              Df) *rgbIPTim[[All, All, 3]];

rgbIPTadapt = rgbIPTim*0;

{rgbIPTadapt[[All, All, 1]] = Rc,  
        rgbIPTadapt[[All, All, 2]] = Gc,  
        rgbIPTadapt[[All, All, 3]] = Bc};

xyzIPTadapt = rgbIPTadapt . Transpose[iCAT02];

dcIPTimage = xyzIPTadapt . Transpose[xyzRGBm];

dcIPTimage = Map[threshold,  dcIPTimage,  {3}];

Max[lmsIPTimage]

<< Statistics`DescriptiveStatistics`

maRed = 
      Quantile[Flatten[dcIPTimage[[All, All, 1]]],  .98]

miRed = Min[dcIPTimage[[All, All, 1]]];

maGreen = 
      Quantile[Flatten[dcIPTimage[[All, All, 2]]],   .98]

miGreen = Min[dcIPTimage[[All, All, 2]]];

maBlue = 
      Quantile[Flatten[dcIPTimage[[All, All, 3]]],   .98]
miBlue = Min[dcIPTimage[[All, All, 3]]];

mightymax  =  Max[{maBlue,  maRed,  maGreen}]

scaleRGB = dcIPTimage * 0 ;

scaleRGB[[All, All, 
            1]] = ((dcIPTimage[[All, All, 1]] - miRed)/(mightymax - miRed));

scaleRGB[[All, All, 
            2]] = ((dcIPTimage[[All, All, 2]] - miGreen)/(mightymax - miGreen));

scaleRGB[[All, All, 
            3]] = ((dcIPTimage[[All, All, 3]] - miBlue)/(mightymax - miBlue));

Show[
        Graphics[
          RasterArray[
            Apply[(RGBColor[#1, #2, #3] &), 
(scaleRGB^((1/1.8))), {2}]]],  AspectRatio -> Automatic];