# extract the exposure gradient for silent pictures 

# read the two raw images
system('dcraw -4 -D regular.CR2 silent.DNG');
regular = double(imread('regular.pgm')) - 1024;
silent = double(imread('silent.pgm')) - 1024;

# for some reason, the CR2 is one line/column larger
regular = regular(1:end-1, 1:end-1);

# extract the green1 channel, to keep things simple
regular = regular(1:2:end, 2:2:end);
silent  = silent (1:2:end, 2:2:end);

# mask out bad pixels
regular(regular < 0) = NaN;
silent(silent < 0) = NaN;

# go to EV space
regular = log2(regular);
silent = log2(silent);

# exposure difference
expo_delta = silent - regular;

# show the images
imshow(regular, []);
set(gca,'position', [0 0 1 1])
print -dpng "-S300,200" regular-ev.png

imshow(silent, []);
set(gca,'position', [0 0 1 1])
print -dpng "-S300,200" silent-ev.png

imshow(expo_delta, []);
set(gca,'position', [0 0 1 1])
print -dpng "-S300,200" expo-delta.png

# plot the exposure for each row, in EV
figure
expo_delta(~isfinite(expo_delta)) = 0;
row_exposure = median(expo_delta');
row = (1:length(row_exposure)) * 2;
plot(row, row_exposure, 'LineWidth', 2)
title('Silent picture gradient')
xlabel('Image row');
ylabel('Exposure difference (EV)')
grid on, axis tight
print -dpng -FHelvetica:14 -r75 row_exposure.png

# plot the exposure for each row, in linear units
# the reference image was taken at 1/4 seconds
figure
exposure_time = 2.^row_exposure * 1/4;
plot(row, exposure_time * 1000, 'LineWidth', 2)
title('Silent picture gradient')
xlabel('Image row');
ylabel('Exposure time (milliseconds)')
grid on, axis([0 max(row) 0 max(exposure_time)*1000])
print -dpng -FHelvetica:14 -r75 row_exposure_time.png

# compute exposure range and rolling shutter
p = polyfit(row(10:end-10), exposure_time(10:end-10), 1);
fit_expo_time = p(2) + p(1) * row;
hold on, plot(row, fit_expo_time * 1000, 'r')

disp(sprintf('Exposure range  : 1/%.2f ... 1/%.2f', 1/fit_expo_time(1), 1/fit_expo_time(end)))
disp(sprintf('Rolling shutter : %.2f seconds', fit_expo_time(end) -fit_expo_time(1)))