Optimal gamma curve for HDMI recording

When converting 12-bit raw data to 8-bit, for transferring it via HDMI, we want to lose as little useful data as possible. Can we find a curve that minimizes the loss?

We already know the data is corrupted by random noise, which has two main components: read noise and shot noise:

noise=read_noise2+shot_noise2

One way to optimize our curve would be to throw away more bits of data in those areas where the noise is higher. If we are lucky, we'll probably throw away only data below the noise floor. If we are not, we'll have to find out what else we can throw away...

Recommended reading:

I will use Octave 4.0 for this experiment.

Noise curve

From the datasheet and experiments, we have found out that read noise, on a dark frame, is about 4 DN, and we know that, at gain x1, a sensel saturates at 13500 electrons.

In [1]:
read_noise = 4;              % 12-bit DN
full_well = 13500;           % electrons

The sensor gains and offsets are set in such a way that black level is around 100 and white level is around 4000. Therefore, the signal range is around 3900 DN.

Armed with this knowledge, we can write the noise model:

In [2]:
function noise = noise_model(dn, read_noise, full_well)
    dn = min(max(dn, 1e-5), 3900);
    gain = full_well / 3900;
    electrons = dn * gain;
    shot_snr = sqrt(electrons);
    shot_noise = dn ./ shot_snr;
    noise = sqrt(read_noise.^2 + shot_noise.^2);
end

In this analysis, the negative values should be interpreted as noise below black level. These values should not be clipped. Therefore, the range of the analyzed digital numbers will be -100:3900. The values would be 0-4000 in the raw data, with a black level of 100 - this roughly match our sensor setup.

We will need a small helper function to set this range on the graphs.

In [3]:
% helper to force the X range on the graphs
function xaxis(xr)
    a = axis;
    a(1:2) = xr;
    axis(a);
end

function xrange12()
    xaxis([-100 3900]);
end

Noise curve:

In [4]:
dn = -100:3900;
noise = noise_model(dn, read_noise, full_well);
subplot(121)
plot(dn, noise)
grid on
set(gca, 'xtick', linspace(0, 4000, 5))
xlabel('Signal (DN)')
ylabel('Noise (DN)')
xrange12()
p = get(gca, 'position'); p(2) = 0.25; p(4) = 0.7; set(gca, 'position', p)
subplot(122)
plot(log2(dn(dn>0)), log2(dn(dn>0)) - log2(noise(dn>0)))
grid on
xlabel('Signal (EV)')
ylabel('SNR (EV)')
p = get(gca, 'position'); p(2) = 0.25; p(4) = 0.7; set(gca, 'position', p)

Okay, so if the noise is about 4 DN in deep shadows, and about 33 DN in highlights, we already know what to start throwing away :)

If the slope (derivative) of our curve will be about 8 times as steep in shadows than in highlights, the quantization error will have the same magnitude relative to noise. Intuition says the ideal curve is the integral of the reciprocal of the noise curve.

In [5]:
lut = cumsum(1./noise);
lut = lut * 255 / lut(end);
plot(dn, lut)
grid on
axis([-100 3900 0 255]),
set(gca, 'ytick', [0 64 128 192 255])

Let's check the quantization error, expressed in noise stdevs.

In [6]:
% reverse curve
rlut = interp1(lut, dn, 0:255);
rlut(1) = dn(1);
plot(0:255, rlut)
grid on
In [7]:
% quantization error: if we alter the 8-bit data by 0.5 (peak roundoff error), how much does the 12-bit data change?
qe = diff(rlut) / 2;
qe(end+1) = qe(end);
plot(rlut, qe ./ interp1(dn,noise,rlut));
axis([-100 3900 0 1])
grid on

So, we managed to get a quantization error of max 0.45 noise stdevs. Not that bad!

Let's verify the result by brute force.

In [8]:
data_8bit = round(lut);
plot(dn, data_8bit)
axis([-100 3900 0 255]),
set(gca, 'ytick', [0 64 128 192 255])
grid on
In [9]:
recovered_12bit = rlut(data_8bit+1);
qe = recovered_12bit - dn;
plot(dn, qe)
grid on
xrange12()

Does the quantization error remain constant relative to the noise levels?

In [10]:
plot(dn, qe ./ noise)
xrange12()
grid on

Answer: yes!

Let's try on a simulated image.

In [11]:
signal = ones(100,1) * (-100:5:3900);
noise = randn(size(signal)) .* noise_model(signal, read_noise, full_well);
image = round(signal + noise);
image = min(max(image + 100, 0), 4095);
imshow(image, [0 4095]);
In [12]:
image_8bit = round(lut(min(max(image + 1, 1), 4001)));
image_recovered12bit = rlut(image_8bit+1);
imshow([image_8bit*16; image_recovered12bit], [0 4095])

Let's see how well we have recovered the original signal.

In [13]:
% difference between original signal and recovered image
err = image_recovered12bit - signal;

% difference relative to (known) noise stdev - should combine the effects of the original noise and the quantization error
errn = err ./ noise_model(signal, read_noise, full_well);

% random noise relative to noise stdev - should be a Gaussian with stdev=1
nn = noise ./ noise_model(signal, read_noise, full_well);

imshow([err; errn*10; nn*10],[])
In [14]:
noise_normalized = std(nn(:))
noise_after_recovery_normalized = std(errn(:))
extra_noise = log2(noise_after_recovery_normalized / noise_normalized)
noise_normalized =  1.0008
noise_after_recovery_normalized =  1.0314
extra_noise =  0.043436

So, going to 8 bits and back introduces only 0.05 stops of extra noise?! Sounds like it's almost lossless!

What the...

The 16-235 trick

There's one more trick to take care of: the HDMI recorder clips everything under 16 and above 235.

In [15]:
function im = scale_hdmi(im)
    im = min(max(im,16),235);
    im = (im - 16) * 255 / (235 - 16);
end

No big deal - instead of our LUT being from 0-255, it will be from 16-235.

In [16]:
noise = noise_model(dn, read_noise, full_well);
lut = cumsum(1./noise);
lut = lut * (235 - 16) / lut(end) + 16;

% we will reuse this, so make it a function
function check_lut(lut, dn, noise)
    rlut = interp1(lut, dn, 16:235);
    rlut(1) = dn(1);

    data_8bit = min(max(round(lut), 16), 235);
    recovered_12bit = rlut(data_8bit-16+1);
    qe = recovered_12bit - dn;
    plot(dn, qe ./ noise)
    xrange12()
    grid on
end

check_lut(lut, dn, noise)

As expected, the roundoff error increased from about 0.45 to about 0.5.

In [17]:
signal = ones(100,1) * (-100:5:3900);
noise = randn(size(signal)) .* noise_model(signal, read_noise, full_well);
image = round(signal + noise);
image = min(max(image + 100, 0), 4095);
err = image_recovered12bit - signal;
errn = err ./ noise_model(signal, read_noise, full_well);
nn = noise ./ noise_model(signal, read_noise, full_well);
noise_normalized = std(nn(:))
noise_after_recovery_normalized = std(errn(:))
extra_noise = log2(noise_after_recovery_normalized / noise_normalized)
noise_normalized =  0.99618
noise_after_recovery_normalized =  1.0314
extra_noise =  0.050086

Right - the roundoff errors introduced about 0.05 stops of extra "noise".

Approximating the LUT with a gamma curve

For practical reasons - that's how the Beta is configured at the moment - we can approximate the optimal LUT curve using 3 parameters:

  • gamma
  • gain
  • offset

It's not really needed, as the FPGA accepts any kind of LUT, not just gamma, but it's easier to configure it that way :P

In [18]:
gamma_lut = @(x, gamma, gain, offset) min(max(((((x + offset) .* gain) / 4095) .^ gamma) * 255, 0), 255);
In [19]:
params = fminsearch(@(p) norm(lut - gamma_lut(dn, p(1), p(2), p(3))), [0.5 1 50])

plot(dn, lut, 'b', dn, gamma_lut(dn, params(1), params(2), params(3)), 'r'), grid on
xrange12()
params =

     0.47189     0.85879   104.01745

We don't really need to represent values below black - 30, so let's optimize the curve for those values only.

In [20]:
dn = -30:3970;
noise = noise_model(dn, read_noise, full_well);
lut = cumsum(1./noise);
lut = lut * (235 - 16) / lut(end) + 16;
In [21]:
params = fminsearch(@(p) norm(lut - gamma_lut(dn, p(1), p(2), p(3))), [0.5 1 50])

plot(dn, lut, 'b', dn, gamma_lut(dn, params(1), params(2), params(3)), 'r'), grid on
xrange12()
params =

    0.51369    0.87105   50.02030

Let's pick some round values:

In [22]:
params = [0.5, 0.9, 50];
plot(dn, lut, 'b', dn, gamma_lut(dn, params(1), params(2), params(3)), 'r'), grid on
xrange12()

Let's check the range of output values:

In [23]:
in = [-20 0 20 500 1000 2000 3000 3800 3900 4000];
round(scale_hdmi(gamma_lut(in, params(1), params(2), params(3))))
ans =

     5    12    18    85   124   181   224   254   255   255

It might clip some highlights, so let's turn the gain down a little bit:

In [24]:
params = [0.5, 0.85, 50];
plot(dn, lut, 'b', dn, gamma_lut(dn, params(1), params(2), params(3)), 'r'), grid on, xrange12()
round(scale_hdmi(gamma_lut(in, params(1), params(2), params(3))))
ans =

     5    12    17    82   120   175   218   247   250   254

Let's check the rounding error once again.

In [25]:
new_lut = gamma_lut(dn, params(1), params(2), params(3));
check_lut(new_lut, dn, noise)

All good!

Higher gains

Can these parameters be used for higher gains?

In [26]:
% very rough values for gain=4
read_noise = 4*3;              % 12-bit DN
full_well = 13500/4;           % electrons
In [27]:
noise_hi = noise_model(dn, read_noise, full_well);
check_lut(new_lut, dn, noise_hi)