# Bits & exposures¶

Experiments based on ideas presented in Grossberg_CPMCV03 (article suggested by @garry23).

## Experiment 1¶

We prove that even without linearizing the camera, simple summation preserves all the information contained in the set of individual exposures. (emphasis added)

Theorem 1:

The sum of a set of images of a scene taken at different exposures includes all the information in the individual exposures.

Test setup:

• linear ramp: 1024 levels between 0 and 1
• 2-bit linear camera
• baseline exposures: 1, 2, 4
• stacking algorithms tested:
1. exposure matching, followed by selecting the non-clipped areas from the brightest input image (easiest)
2. exposure matching, followed by selecting the brightest pixel out of all of the stacked images (i.e. what I'm doing in ceronoice)
3. exposure matching, followed by averaging of non-clipped areas (i.e. what I wanted to implement in ceronoice)
4. simple summation (as suggested in the paper), followed by straightening the response curve
In :
%plot -s 2000,100

# create a linear ramp using 1024 values between 0 and 1
scene = ones(50, 1) * linspace(0, 1, 1024);

# show the ramp image (test scene)
figure, imshow(scene);

# some displays may not display a smooth ramp, but a posterized one
# adding a bit of noise solves the issue
figure, imshow(scene + randn(size(scene))/200)

# displaying with gamma correction will reveal shadow details
figure, imshow(sqrt(scene))

# noise to avoid posterization
figure, imshow(sqrt(max(scene + randn(size(scene))/200, 0)))    In :
function im = camera(scene, bits, exposure)
# multiply scene by linear exposure
# clip everything below 0 (black) and 1 (white)
input = min(max(scene * exposure, 0), 1);

# simulate a n-bit camera
levels = 2 ^ bits;
im = round(input * levels) / levels;
end

In :
%plot -s 950,950

# simulate bracketing using 3 "baseline" exposures and a 2-bit camera
im1 = camera(scene, 2, 1);
im2 = camera(scene, 2, 2);
im4 = camera(scene, 2, 4);

# combine the 3 exposures
# and replace overexposed areas from darker images
rep = im4 / 4;
rep(im4 == 1) = im2(im4 == 1) / 2;
rep(im2 == 1) = im1(im2 == 1);

# another algorithm: just take the max value from all images after exposure matching
# this is what I'm doing in ceronoice
crn = max(max(im1, im2/2), im4/4);

# another algorithm: average all the exposures covering each tonal range,
# after matching their brightness (i.e. scaling by 1 / exposure)
av1 = im1;
av1(im2 < 1) = (im1(im2 < 1) + im2(im2 < 1) / 2) / 2;
av1(im4 < 1) = (im1(im4 < 1) + im2(im4 < 1) / 2 + im4(im4 < 1) / 4) / 3;

# plain average (i.e. simple summation followed by normalization)
# as we are working on floating point, the division by 3 won't make any difference in bit depth
avg = (im1 + im2 + im4) / 3;

# simple summation of the exposure, followed by piecewise straightening of the response curve

# find the response curve of the plain average (yes, I'm too lazy to do the math)
ramp = linspace(0, 1, 100000);
ravg = (ramp + min(2*ramp, 1) + min(4*ramp, 1)) / 3;

# undo that response curve
av2 = interp1(ravg, ramp, avg);

figure, imshow(sqrt([ im1;
im2;
im4;
rep;
crn;
av1;
avg;
av2;
scene]))

figure, hold on
plot(mean(scene), mean(im1), 'k:'),
plot(mean(scene), mean(im2)/2, 'b:'),
plot(mean(scene), mean(im4)/4, 'r--'),
plot(mean(scene), mean(rep), 'k.-')
plot(mean(scene), mean(crn), 'y.-')
plot(mean(scene), mean(av1), 'm.-')
plot(mean(scene), mean(avg), 'b.-'),
plot(ramp, ravg, 'b--')
plot(mean(scene), mean(av2), 'g.-')
plot(mean(scene), mean(scene), 'r')

# re-plot these to make them a little more visible (i.e. bring to front)
plot(mean(scene), mean(rep), 'k--')
plot(mean(scene), mean(im1), 'k:'),
plot(mean(scene), mean(im2)/2, 'b:'),
plot(mean(scene), mean(im4)/4, 'r--'),

legend(...
'Single exposure (1.0)', ...
'Single exposure (2.0)', ...
'Single exposure (4.0)', ...
'Replacement of overexposed areas', ...
'Max of brightness-matched images (ceronoice)', ...
'Average of non-clipped areas after matching brightness', ...
'Plain average: (im1 + im2 + im4) / 3', ...
'Response curve of plain average', ...
'Corrected average (simple summation + piecewise straightening)', ...
'Ground truth', ...
'location', 'southeast')

title(sprintf([ ...
'Image stacking algorithms at low bit depths\n\n' ...
'2-bit linear camera' ]))
axis square, axis equal, axis([0 1 0 1])  Results:

1. simple summation (as suggested in the paper), followed by straightening the response curve (best)
2. exposure matching, followed by averaging of non-clipped areas (not bad, but has some issues with nonuniform spacing of levels in deep shadows)
3. exposure matching, followed by selecting the non-clipped areas from the brightest input image (it picks the best bits from each exposure, but doesn't improve upon that)
4. exposure matching, followed by selecting the brightest pixel out of all of the stacked images (output biased towards brighter levels, otherwise very similar to 3; there may be differences with noisy input images)

TLDR: it is important to perform exposure corrections after a simple summation of the input images.

TLDR: the paper authors were right. So easy to overlook such as simple issue.

## Experiment 2¶

Adding Bits to a Linear Camera

We start with the simple example of combining 3 exposures from a linear real camera with 4 bits. Our algorithm specifies exposures close together (1, 1.05, 1.11). The response from the resulting effective camera is shown with the response resulting from the baseline exposures (1, 2, 4) in Fig. 5(a). Using the exposures computed with our algorithm allows us to achieve a more uniform response.

Test setup:

• same linear ramp (1024 levels between 0 and 1)
• 4-bit camera (16 levels captured)
• compare the two exposure sequences (1, 1.05, 1.11) and (1, 2, 4) after stacking them with the best algorithm from previous test.
In :
%plot -s 950,950

# "optimal" exposures (according to paper authors)
ima = camera(scene, 4, 1);
imb = camera(scene, 4, 1.05);
imc = camera(scene, 4, 1.11);

# "baseline" exposures
im1 = camera(scene, 4, 1);
im2 = camera(scene, 4, 2);
im4 = camera(scene, 4, 4);

# stack the "optimal" exposures by simple summation followed by curve straightening
avg_o = (ima + imb + imc) / 3;
ramp_o = linspace(0, 1, 100000);
ravg_o = (ramp_o + min(ramp_o * 1.05, 1) + min(ramp_o * 1.11, 1)) / 3;
av2_o = interp1(ravg_o, ramp_o, avg_o);

# stack the "baseline" exposures in the same way
avg_b = (im1 + im2 + im4) / 3;
ramp_b = linspace(0, 1, 100000);
ravg_b = (ramp_b + min(ramp_b * 2, 1) + min(ramp_b * 4, 1)) / 3;
av2_b = interp1(ravg_b, ramp_b, avg_b);

figure
imshow(sqrt([ ima;
imb;
imc;
avg_o;
avg_b;
av2_o;
av2_b;
scene ]))

figure, hold on
plot(mean(scene), mean(ima), 'k'),
plot(mean(scene), mean(avg_o), 'b'),
plot(mean(scene), mean(avg_b), 'm'),
plot(ramp_o, ravg_o, 'b--')
plot(ramp_b, ravg_b, 'm--')
plot(mean(scene), mean(av2_o), 'g.-')
plot(mean(scene), mean(av2_b), 'r.-')
plot(mean(scene), mean(scene), 'r', 'linewidth', 2)

# re-plot the single exposure to bring it to front
plot(mean(scene), mean(ima), 'k'),

legend(...
'Single exposure (1.0)', ...
'Plain average of optimal exposures: (ima + imb + imc) / 3', ...
'Plain average of baseline exposures: (im1 + im2 + im4) / 3', ...
'Response curve of plain average of optimal exposures', ...
'Response curve of plain average of baseline exposures', ...
'Corrected average of optimal exposures', ...
'Corrected average of baseline exposures', ...
'Ground truth', ...
'location', 'southeast')

title(sprintf([ ...
'Adding bits to a linear camera\n\n' ...
'4-bit linear camera (simulated)\n' ...
'Optimal exposures: 1, 1.05, 1.11\n' ...
'Baseline exposures: 1, 2, 4' ]))
axis square, axis equal, axis([0 1 0 1])  Results:

• baseline exposures (1, 2, 4): bit depth much better in shadows, no improvement in highlights when compared to a single exposure (red curve)
• optimal exposures (1, 1.05, 1.11): bit depth significantly better in highlights, but shadows are not much better, compared to a single exposure (green curve)

The suggested exposure sequence (1, 1.05, 1.11) does add some "bits" in midtones and highlights, compared to the baseline sequence (1, 2, 4).

Exercise for the reader: count the extra "bits" :)