This example shows how to estimate illumination and perform white balance of a scene using three different illumination algorithms. The example compares the estimated illuminant to a ground truth illuminant calculated using an X-Rite® ColorChecker® chart.

Eyes are very good at judging what is white under different lighting conditions. Digital cameras, however, without some kind of adjustment, can easily capture unrealistic images with a strong color cast. Automatic white balance (AWB) algorithms try to correct for the ambient light with minimum input from the user, so that the resulting image looks like what our eyes would see.

Automatic white balancing is done in two steps:

Step 1: Estimate the scene illuminant.

Step 2: Correct the color balance of the image.

Several different algorithms exist to estimate scene illuminant. The performance of each algorithm depends on the scene, lighting, and imaging conditions. This example judges the quality of three algorithms for illuminant estimation for one specific image by comparing them to the ground truth scene illuminant:

After the ambient light is known, then correcting the colors in the image (Step 2) is an easy and fixed process.

AWB algorithms are usually applied on the raw image data after a minimal amount of preprocessing, before the image is compressed and saved to the memory card.

Read a 16-bit raw image into the workspace. `foosballraw.tiff`

is an image file that contains raw sensor data after correcting the black level and scaling the intensities to 16 bits per pixel. This image is free of the white balancing done by the camera, as well as other preprocessing operations such as demosaicing, denoising, chromatic aberration compensation, tone adjustments, and gamma correction.

`A = imread('foosballraw.tiff');`

Digital cameras use a color filter array superimposed on the imaging sensor to simulate color vision, so that each pixel is sensitive to either red, green or blue. To recover the missing color information at every pixel, interpolate using the `demosaic`

function. The Bayer pattern used by the camera with which the photo was captured (Canon EOS 30D) is RGGB.

`A = demosaic(A,'rggb');`

The image `A`

contains linear RGB values. Linear RGB values are appropriate for estimating scene illuminant and correcting the color balance of an image. However, if you try to display the linear RGB image, it will appear very dim, because of the nonlinear characteristic of display devices. Therefore, for display purposes, gamma-correct the image to the sRGB color space using the `lin2rgb`

function.

A_sRGB = lin2rgb(A);

Display the demosiaced image before and after gamma correction.

```
montage({A,A_sRGB})
title('Original Image Before and After Gamma Correction')
```

Calculate the ground truth illuminant using the X-Rite ColorChecker chart that is included in the scene. This chart consists of 24 neutral and color patches with known spectral reflectances.

Detect the chart in the gamma-corrected image by using the `colorChecker`

function. The linear RGB image is too dark for `colorChecker`

to detect the chart automatically.

chart_sRGB = colorChecker(A_sRGB);

Confirm that the chart is detected correctly.

displayChart(chart_sRGB)

Get the coordinates of the registration points at the four corners of the chart.

registrationPoints = chart_sRGB.RegistrationPoints;

Create a new `colorChecker`

object from the linear RGB data. Specify the location of the chart using the coordinates of the registration points.

`chart = colorChecker(A,"RegistrationPoints",registrationPoints);`

Measure the RGB values of all color patches using the `measureColor`

function.

colors = measureColor(chart)

`colors=`*24×9 table*
ROI Color Measured_R Measured_G Measured_B Reference_L Reference_a Reference_b Delta_E
___ ________________ __________ __________ __________ ___________ ___________ ___________ _______
1 {'DarkSkin' } 1773 2284 1106 37.54 14.37 14.92 40.763
2 {'LightSkin' } 5691 7627 4245 64.66 19.27 17.5 61.046
3 {'BlueSky' } 1919 5214 4677 49.32 -3.82 -22.54 49.222
4 {'Foliage' } 1576 3376 1216 43.46 -12.74 22.72 46.172
5 {'BlueFlower' } 2997 6253 5997 54.94 9.61 -24.79 55.312
6 {'BluishGreen' } 3676 12020 7388 70.48 -32.26 -0.37 56.895
7 {'Orange' } 6250 5276 1213 62.73 35.83 56.5 82.432
8 {'PurplishBlue'} 1226 3779 5302 39.43 10.75 -45.17 55.714
9 {'ModerateRed' } 4581 3277 2074 50.57 48.64 16.67 67.971
10 {'Purple' } 1233 1767 1889 30.1 22.54 -20.87 41.882
11 {'YellowGreen' } 4939 10555 2657 71.77 -24.13 58.19 72.055
12 {'OrangeYellow'} 7410 8522 1608 71.51 18.24 67.37 83.134
13 {'Blue' } 567 2300 3888 28.37 15.42 -49.8 55.878
14 {'Green' } 2087 6486 2278 54.38 -39.72 32.27 61.962
15 {'Red' } 3702 1986 968 42.43 51.05 28.62 68.861
16 {'Yellow' } 8906 12782 2422 81.8 2.67 80.41 87.375
⋮

To measure the ground truth illuminant, only use the 6 neutral patches on the bottom row of the chart. Extract the RGB values of the six neutral (achromatic) patches on the bottom row.

grayPatchRGB = colors{19:end,{'Measured_R','Measured_G','Measured_B'}}; grayPatchRGB = im2double(grayPatchRGB);

The ground truth illuminant is the average color of the neutral patches.

illuminant_groundtruth = mean(grayPatchRGB)

`illuminant_groundtruth = `*1×3*
0.0693 0.1423 0.0943

When testing the AWB algorithms, prevent the algorithms from unfairly taking advantage of the chart by masking out the chart.

Create a polygon ROI over the chart by using the `drawpolygon`

function. Specify the vertices of the polygon as the registration points.

`chartROI = drawpolygon("Position",registrationPoints);`

Convert the polygon ROI to a binary mask by using the `createMask`

function.

mask_chart = createMask(chartROI);

Invert the mask. Pixels within the chart are excluded from the mask and pixels of the rest of the scene are included in the mask.

mask_scene = ~mask_chart;

To confirm the accuracy of the mask, display the mask over the image. Pixels included in the mask have a blue tint.

imshow(labeloverlay(A_sRGB,mask_scene));

You can consider an illuminant as a vector in 3-D RGB color space. The magnitude of the estimated illuminant does not matter as much as its direction, because the direction of the illuminant is what is used to white balance an image.

To evaluate the quality of an estimated illuminant, compute the angular error between the estimated illuminant and the ground truth. Angular error is the angle (in degrees) formed by the two vectors. The smaller the angular error, the better the estimation is.

To better understand the concept of angular error, consider the following visualization of an arbitrary illuminant and the ground truth measured using the ColorChecker chart. The `plotColorAngle`

helper function plots a unit vector of an illuminant in 3-D RGB color space, and is defined at the end of the example.

sample_illuminant = [0.066 0.1262 0.0691]; p = plot3([0 1],[0 1],[0,1],'LineStyle',':','Color','k'); ax = p.Parent; hold on plotColorAngle(illuminant_groundtruth,ax) plotColorAngle(sample_illuminant,ax) title('Illuminants in RGB space') view(28,36) legend('Achromatic Line','Ground Truth Illuminant','Sample Illuminant') grid on axis equal

The White Patch Retinex algorithm for illuminant estimation assumes that the scene contains a bright achromatic patch. This patch reflects the maximum light possible for each color band, which is the color of the scene illuminant. Use the `illumwhite`

function to estimate illumination using the White Patch Retinex algorithm.

Estimate the illuminant using all the pixels in the scene. Exclude the ColorChecker chart from the scene by using the `'Mask'`

name-value pair argument.

```
percentileToExclude = 0;
illuminant_wp1 = illumwhite(A,percentileToExclude,'Mask',mask_scene);
```

Compute the angular error for the illuminant estimated with White Patch Retinex.

```
err_wp1 = colorangle(illuminant_wp1, illuminant_groundtruth);
disp(['Angular error for White Patch with percentile=0: ' num2str(err_wp1)])
```

Angular error for White Patch with percentile=0: 16.5377

White balance the image using the `chromadapt`

function. Specify the estimated illuminant and indicate that color values are in the linear RGB color space.

B_wp1 = chromadapt(A,illuminant_wp1,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image.

```
B_wp1_sRGB = lin2rgb(B_wp1);
figure
imshow(B_wp1_sRGB)
title('White-Balanced Image using White Patch Retinex with percentile=0')
```

The White Patch Retinex algorithm does not perform well when pixels are overexposed. To improve the performance of the algorithm, exclude the top 1% of the brightest pixels.

```
percentileToExclude = 1;
illuminant_wp2 = illumwhite(A,percentileToExclude,'Mask',mask_scene);
```

Calculate the angular error for the estimated illuminant. The error is less than when estimating the illuminant using all pixels.

```
err_wp2 = colorangle(illuminant_wp2,illuminant_groundtruth);
disp(['Angular error for White Patch with percentile=1: ' num2str(err_wp2)])
```

Angular error for White Patch with percentile=1: 5.032

White balance the image in the linear RGB color space using the estimated illuminant.

B_wp2 = chromadapt(A,illuminant_wp2,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image with the new illuminant.

```
B_wp2_sRGB = lin2rgb(B_wp2);
imshow(B_wp2_sRGB)
title('White-Balanced Image using White Patch Retinex with percentile=1')
```

The Gray World algorithm for illuminant estimation assumes that the average color of the world is gray, or achromatic. Therefore, it calculates the scene illuminant as the *average* RGB value in the image. Use the `illumgray`

function to estimate illumination using the Gray World algorithm.

First, estimate the scene illuminant using all pixels of the image, excluding those corresponding to the ColorChecker chart. The `illumgray`

function provides a parameter to specify the percentiles of bottom and top values (ordered by brightness) to exclude. Here, specify the percentiles as 0.

```
percentileToExclude = 0;
illuminant_gw1 = illumgray(A,percentileToExclude,'Mask',mask_scene);
```

Calculate the angular error between the estimated illuminant and the ground truth illuminant.

```
err_gw1 = colorangle(illuminant_gw1,illuminant_groundtruth);
disp(['Angular error for Gray World with percentiles=[0 0]: ' num2str(err_gw1)])
```

Angular error for Gray World with percentiles=[0 0]: 5.0414

White balance the image in the linear RGB color space using the estimated illuminant.

B_gw1 = chromadapt(A,illuminant_gw1,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image.

```
B_gw1_sRGB = lin2rgb(B_gw1);
imshow(B_gw1_sRGB)
title('White-Balanced Image using Gray World with percentiles=[0 0]')
```

The Gray World algorithm does not perform well when pixels are underexposed or overexposed. To improve the performance of the algorithm, exclude the top 1% of the darkest and brightest pixels.

```
percentileToExclude = 1;
illuminant_gw2 = illumgray(A,percentileToExclude,'Mask',mask_scene);
```

Calculate the angular error for the estimated illuminant. The error is less than when estimating the illuminant using all pixels.

```
err_gw2 = colorangle(illuminant_gw2, illuminant_groundtruth);
disp(['Angular error for Gray World with percentiles=[1 1]: ' num2str(err_gw2)])
```

Angular error for Gray World with percentiles=[1 1]: 5.1092

White balance the image in the linear RGB color space using the estimated illuminant.

B_gw2 = chromadapt(A,illuminant_gw2,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image with the new illuminant.

```
B_gw2_sRGB = lin2rgb(B_gw2);
imshow(B_gw2_sRGB)
title('White-Balanced Image using Gray World with percentiles=[1 1]')
```

Cheng's illuminant estimation method draws inspiration from spatial domain methods such as Grey Edge [4], which assumes that the gradients of an image are achromatic. They show that Grey Edge can be improved by artificially introducing strong gradients by shuffling image blocks, and conclude that the strongest gradients follow the direction of the illuminant. Their method consists in ordering pixels according to the norm of their projection along the direction of the mean image color, and retaining the bottom and top percentile. These two groups correspond to strong gradients in the image. Finally, they perform a *principal component analysis* (PCA) on the retained pixels and return the first component as the estimated illuminant. Use the `illumpca`

function to estimate illumination using Cheng's PCA algorithm.

First, estimate the illuminant using the default percentage value of Cheng's PCA method, excluding those corresponding to the ColorChecker chart.

`illuminant_ch2 = illumpca(A,'Mask',mask_scene);`

Calculate the angular error between the estimated illuminant and the ground truth illuminant.

```
err_ch2 = colorangle(illuminant_ch2,illuminant_groundtruth);
disp(['Angular error for Cheng with percentage=3.5: ' num2str(err_ch2)])
```

Angular error for Cheng with percentage=3.5: 5.0159

White balance the image in the linear RGB color space using the estimated illuminant.

B_ch2 = chromadapt(A,illuminant_ch2,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image.

```
B_ch2_sRGB = lin2rgb(B_ch2);
imshow(B_ch2_sRGB)
title('White-Balanced Image using Cheng with percentile=3.5')
```

Now, estimate the scene illuminant using the bottom and top 5% of pixels along the direction of the mean color. The second argument of the `illumpca`

function specifies the percentiles of bottom and top values (ordered by brightness) to exclude.

`illuminant_ch1 = illumpca(A,5,'Mask',mask_scene);`

Calculate the angular error between the estimated illuminant and the ground truth illuminant. The error is less than when estimating the illuminant using the default percentage.

```
err_ch1 = colorangle(illuminant_ch1, illuminant_groundtruth);
disp(['Angular error for Cheng with percentage=5: ' num2str(err_ch1)])
```

Angular error for Cheng with percentage=5: 4.7451

White balance the image in the linear RGB color space using the estimated illuminant.

B_ch1 = chromadapt(A,illuminant_ch1,'ColorSpace','linear-rgb');

Display the gamma-corrected white-balanced image.

```
B_ch1_sRGB = lin2rgb(B_ch1);
imshow(B_ch1_sRGB)
title('White-Balanced Image using Cheng with percentage=5')
```

To find the best parameter to use for each algorithm, you can sweep through a range and calculate the angular error for each of them. The parameters of the three algorithms have different meanings, but the similar ranges of these parameters makes it easy to programmatically search for the best one for each algorithm.

param_range = 0:0.25:5; err = zeros(numel(param_range),3); for k = 1:numel(param_range) % White Patch illuminant_wp = illumwhite(A,param_range(k),'Mask',mask_scene); err(k,1) = colorangle(illuminant_wp,illuminant_groundtruth); % Gray World illuminant_gw = illumgray(A,param_range(k),'Mask',mask_scene); err(k,2) = colorangle(illuminant_gw,illuminant_groundtruth); % Cheng if (param_range(k) ~= 0) illuminant_ch = illumpca(A,param_range(k),'Mask',mask_scene); err(k,3) = colorangle(illuminant_ch,illuminant_groundtruth); else % Cheng's algorithm is undefined for percentage=0. err(k,3) = NaN; end end

Display a heatmap of the angular error using the `heatmap`

function. Dark blue colors indicate a low angular error while yellow colors indicate a high angular error. The optimal parameter has the smallest angular error.

heatmap(err,'Title','Angular Error','Colormap',parula(length(param_range)), ... 'XData',["White Patch" "Gray World" "Cheng's PCA"], ... 'YLabel','Parameter Value','YData',string(param_range));

Find the best parameter for each algorithm.

[~,idx_best] = min(err); best_param_wp = param_range(idx_best(1)); best_param_gw = param_range(idx_best(2)); best_param_ch = param_range(idx_best(3)); fprintf('The best parameter for White Patch is %1.2f with angular error %1.2f degrees\n', ... best_param_wp,err(idx_best(1),1));

The best parameter for White Patch is 0.25 with angular error 3.35 degrees

fprintf('The best parameter for Gray World is %1.2f with angular error %1.2f degrees\n', ... best_param_gw,err(idx_best(2),2));

The best parameter for Gray World is 0.00 with angular error 5.04 degrees

fprintf('The best parameter for Cheng is %1.2f with angular error %1.2f degrees\n', ... best_param_ch,err(idx_best(3),3));

The best parameter for Cheng is 0.50 with angular error 1.74 degrees

Calculate the estimated illuminant for each algorithm using the best parameter.

best_illum_wp = illumwhite(A,best_param_wp,'Mask',mask_scene); best_illum_gw = illumgray(A,best_param_gw,'Mask',mask_scene); best_illum_ch = illumpca(A,best_param_ch,'Mask',mask_scene);

Display the angular error of each best illuminant in the RGB color space.

p = plot3([0 1],[0 1],[0,1],'LineStyle',':','Color','k'); ax = p.Parent; hold on plotColorAngle(illuminant_groundtruth,ax) plotColorAngle(best_illum_wp,ax) plotColorAngle(best_illum_gw,ax) plotColorAngle(best_illum_ch,ax) title('Best Illuminants in RGB space') view(28,36) legend('Achromatic Line','Ground Truth','White Patch','Gray World','Cheng') grid on axis equal

Calculate the optimal white-balanced images for each algorithm using the best illuminant.

B_wp_best = chromadapt(A,best_illum_wp,'ColorSpace','linear-rgb'); B_wp_best_sRGB = lin2rgb(B_wp_best); B_gw_best = chromadapt(A,best_illum_gw,'ColorSpace','linear-rgb'); B_gw_best_sRGB = lin2rgb(B_gw_best); B_ch_best = chromadapt(A,best_illum_ch,'ColorSpace','linear-rgb'); B_ch_best_sRGB = lin2rgb(B_ch_best);

Display the optimal white-balanced images for each algorithm in a montage.

figure montage({B_wp_best_sRGB,B_gw_best_sRGB,B_ch_best_sRGB},'Size',[1 3]) title('Montage of Best White-Balanced Images: White Point, Gray World, Cheng')

This quick shootout between two classic illuminant estimation algorithms and a more recent one shows that Cheng's method, using the top and bottom 0.75% darkest and brightest pixels, wins for that particular image. However, this result should be taken with a grain of salt.

First, the ground truth illuminant was measured using a ColorChecker chart and is sensitive to shot and sensor noise. The ground truth illuminant of a scene can be better estimated using a spectrophotometer.

Second, the ground truth illuminant is estimated as the mean color of the neutral patches. It is common to use the median instead of the mean, which could shift the ground truth by a significant amount. For example, for the image in this study, using the same pixels, the median color and the mean color of the neutral patches are 0.5 degrees apart, which in some cases can be more than the angular error of the illuminants estimated by different algorithms.

Third, a full comparison of illuminant estimation algorithms should use a variety of images taken under different conditions. One algorithm might work better than the others for a particular image, but might perform poorly over the entire data set.

The `plotColorAngle`

function plots a unit vector of an illuminant in 3-D RGB color space. The input argument `illum`

specifies the illuminant as an RGB color and the input argument `ax`

specifies the axes on which to plot the unit vector.

function plotColorAngle(illum,ax) R = illum(1); G = illum(2); B = illum(3); magRGB = norm(illum); plot3([0 R/magRGB],[0 G/magRGB],[0 B/magRGB], ... 'Marker','.','MarkerSize',10,'Parent',ax) xlabel('R') ylabel('G') zlabel('B') xlim([0 1]) ylim([0 1]) zlim([0 1]) end

[1] Ebner, Marc. White Patch Retinex, Color Constancy. John Wiley & Sons, 2007. ISBN 978-0-470-05829-9.

[2] Ebner, Marc. The Gray World Assumption, Color Constancy. John Wiley & Sons, 2007. ISBN 978-0-470-05829-9.

[3] Cheng, Dongliang, Dilip K. Prasad, and Michael S. Brown. "Illuminant estimation for color constancy: why spatial-domain methods work and the role of the color distribution." JOSA A 31.5 (2014): 1049-1058.

[4] Van De Weijer, Joost, Theo Gevers, and Arjan Gijsenij. "Edge-based color constancy." IEEE Transactions on image processing 16.9 (2007): 2207-2214.

`chromadapt`

| `colorangle`

| `colorChecker`

| `illumgray`

| `illumpca`

| `illumwhite`

| `lin2rgb`

| `measureColor`

| `rgb2lin`