Flat frame experiment¶

Long version - for very dark spots in H.264 clips¶

Noticed some artifacts in a few H.264 clips recorded with 5D Mark III and 50/1.8 at f/22, crop_rec set to 1:1. Apparently these are caused by dust on the sensor and made obvious at narrow apertures. Here's one frame:

In [1]:
%%shell
ffmpeg -loglevel quiet -hide_banner -i sample.mov -vframes 1 -y sample.ppm
convert sample.ppm sample.jpg # for preview

sample

Luckily, the spot was large enough to be noticed on the camera screen, so I've recorded a short clip with just the blank sky, moving the camera around to average out the clouds.

In [2]:
%%shell
ffmpeg -loglevel quiet -hide_banner -i sky.mov -f image2pipe -vcodec ppm - | convert - -evaluate-sequence mean sky.ppm
convert sky.ppm sky.jpg # for preview

sky

Let's get a flat frame out of this!

In [3]:
sky = double(imread('sky.ppm'));
sr = sky(:,:,1);
sg = sky(:,:,2);
sb = sky(:,:,3);
ref = sky * 0;
ref(:,:,1) = mean(sr(:));
ref(:,:,2) = mean(sg(:));
ref(:,:,3) = mean(sb(:));
flat = sky ./ ref;
imshow(flat / 2)

Better just make it grayscale.

In [4]:
imshow([flat(:,:,1), flat(:,:,2), flat(:,:,3)])

Keeping just the green channel, as it seems to be cleaner.

In [5]:
gflat = flat(:,:,2);
imshow(gflat)
mean(gflat(:))

ans =  1.0000

LGTM, now let's correct our test frame.

In [6]:
im = double(imread('sample.ppm'));
ca(:,:,1) = im(:,:,1) ./ gflat;
ca(:,:,2) = im(:,:,2) ./ gflat;
ca(:,:,3) = im(:,:,3) ./ gflat;
imshow(ca/255)

LOL, now that's some overcorrection :)

Let's try some gamma correction (trial and error).

In [7]:
g = 1.6;
cb(:,:,1) = (im(:,:,1).^g ./ gflat).^(1/g);
cb(:,:,2) = (im(:,:,2).^g ./ gflat).^(1/g);
cb(:,:,3) = (im(:,:,3).^g ./ gflat).^(1/g);
imshow(cb/255)

Individual channels:

In [8]:
imshow([cb(:,:,1); cb(:,:,2); cb(:,:,3)])

Not perfect, but not that bad either.

Maybe using per-channel corrections could help?

In [9]:
g = 1.6;
cc(:,:,1) = (im(:,:,1).^g ./ flat(:,:,1)).^(1/g);
cc(:,:,2) = (im(:,:,2).^g ./ flat(:,:,2)).^(1/g);
cc(:,:,3) = (im(:,:,3).^g ./ flat(:,:,3)).^(1/g);
imshow(cc/255)
In [10]:
imshow([cc(:,:,1); cc(:,:,2); cc(:,:,3)])

Marginal improvement.

It should work a lot better if you try to correct some raw footage, as the response there is linear. In H.264, it's not. Reversing the curves for the picture style used in this clip might also help.

Let's try something easier. That spot appears on the sky, so we already know the input image, the correction data and we can easily guess the expected result. Zooming in:

In [11]:
%plot -s 900,300
x = im(415:465,800:850,:);
imshow([x(:,:,1), x(:,:,2), x(:,:,3)] / 255)
In [12]:
%plot -s 900,300
f = flat(415:465,800:850,:);
imshow([f(:,:,1), f(:,:,2), f(:,:,3)] / 2)
In [13]:
%plot -s 900,300
r = x;
r(:,:,1) = median(r(:,:,1)(:));
r(:,:,2) = median(r(:,:,2)(:));
r(:,:,3) = median(r(:,:,3)(:));
imshow([r(:,:,1), r(:,:,2), r(:,:,3)] / 255)
In [14]:
plot(f(:,:,1)(:), (x(:,:,1) ./ r(:,:,1))(:), '.r'); hold on
plot(f(:,:,2)(:), (x(:,:,2) ./ r(:,:,2))(:), '.g')
plot(f(:,:,3)(:), (x(:,:,3) ./ r(:,:,3))(:), '.b')

t = linspace(0.2, 1.2, 100);
p = polyfit(f(:), x(:) ./ r(:), 3);
plot(t, polyval(p, t), 'b')

Looks like we might have identified a better correction curve.

In [15]:
%plot -s 900,300
xd(:,:,1) = x(:,:,1) ./ polyval(p, f(:,:,1));
xd(:,:,2) = x(:,:,2) ./ polyval(p, f(:,:,2));
xd(:,:,3) = x(:,:,3) ./ polyval(p, f(:,:,3));
imshow([xd(:,:,1), xd(:,:,2), xd(:,:,3)] / 255)
In [16]:
plot(f(:,:,1)(:), xd(:,:,1)(:), '.r'); hold on
plot(f(:,:,2)(:), xd(:,:,2)(:), '.g')
plot(f(:,:,3)(:), xd(:,:,3)(:), '.b')

Meaning: we probably can't do much better than that.

In [17]:
cd(:,:,1) = im(:,:,1) ./ polyval(p, flat(:,:,1));
cd(:,:,2) = im(:,:,2) ./ polyval(p, flat(:,:,2));
cd(:,:,3) = im(:,:,3) ./ polyval(p, flat(:,:,3));
imshow(cd / 255)

Cheating a bit: let's average the sample movie.

In [18]:
%%shell
ffmpeg -loglevel quiet -hide_banner -i sample.mov -f image2pipe -vcodec ppm - | convert - -evaluate-sequence mean avg.ppm
convert avg.ppm avg.jpg # for preview

avg

We could just take the correction for that troublesome spot from here, as it didn't touch the airplane.

In [19]:
%plot -s 900,300
avg = double(imread('avg.ppm'));
avg = avg(415:465,800:850,:);
aref = avg;
aref(:,:,1) = median(avg(:,:,1)(:));
aref(:,:,2) = median(avg(:,:,2)(:));
aref(:,:,3) = median(avg(:,:,3)(:));
favg = avg ./ aref;
fa = favg .^ g;
h = f; h(1:end/2,:,:) = fa(1:end/2,:,:);
imshow([h(:,:,1), h(:,:,2), h(:,:,3)] / 2)
In [20]:
xe(:,:,1) = (x(:,:,1).^g ./ fa(:,:,1)).^(1/g);
xe(:,:,2) = (x(:,:,2).^g ./ fa(:,:,2)).^(1/g);
xe(:,:,3) = (x(:,:,3).^g ./ fa(:,:,3)).^(1/g);
imshow([xe(:,:,1), xe(:,:,2), xe(:,:,3)] / 255)

Let's replace that spot, taking care not to introduce a square around it.

In [21]:
flat(415:465,800:850,1) = fa(:,:,1) * median(flat(415:465,800:850,1)(:)) / median(fa(:,:,1)(:));
flat(415:465,800:850,2) = fa(:,:,2) * median(flat(415:465,800:850,2)(:)) / median(fa(:,:,2)(:));
flat(415:465,800:850,3) = fa(:,:,3) * median(flat(415:465,800:850,3)(:)) / median(fa(:,:,3)(:));
fchk = flat(405:475, 790:860, :);
imshow([fchk(:,:,1), fchk(:,:,2), fchk(:,:,3)] / 2)
In [22]:
ce(:,:,1) = (im(:,:,1).^g ./ flat(:,:,1)).^(1/g);
ce(:,:,2) = (im(:,:,2).^g ./ flat(:,:,2)).^(1/g);
ce(:,:,3) = (im(:,:,3).^g ./ flat(:,:,3)).^(1/g);
imshow(ce/255)

Happy pixel peeping :)

Correcting the entire footage is left as an exercise for the reader. Need some hints?

OK, now gotta clean that sensor.