% Version tracking: Curated version of DynRgRestoreLG15.m

% This function processes multiple-exposure movies to recover their dynamic
% range approximately back to one of single-exposure movies. The raw movies
% should have linear encoding (gamma=1 in PCC) and their dynamic range
% pre-expanded as much as possible (larger gain setting in PCC). The
% processed movie is automatically saved in avi format with linear
% encoding, grayscale, and uncompressed.

% INPUT variables
% RawName: file name of the raw movie
% ProcName: file name of the processed movie
% Exposures: number of light pulses used for each image
% HeadRoom: interval at the top of the intensity range, reserved to prevent
% saturation of raw pixels brighter than the average background

% OUTPUT variables (besides generating the enhanced movie)
% CombinedDropPattern: average of whole movie (shows the mean drop location)

% EXAMPLE command
% [CombinedDropPattern] = MultiexposureContrastRecovery('Mov025Linear_PARTIAL.avi','M025DropSplit74mm_PARTIAL.avi',12,6);

function [CombinedDropPattern] = MultiexposureContrastRecovery(RawName,ProcName,Exposures,HeadRoom)

% Definitions
LowCut = 0.96;  % To exclude pixels from randomly-located objects (e.g. drops) from movie averages
HiCut = 1.04;   % To exclude pixels from randomly-located objects (e.g. drops) from movie averages

% Start timer
tic

% Load linearly encoded movie file
RawMov = VideoReader(RawName);
vidFrames = read(RawMov);
numFrames = get(RawMov, 'NumberOfFrames');
ColormapBW = [linspace(0,1,256)' linspace(0,1,256)' linspace(0,1,256)'];
fprintf('Movie loaded. '); toc


% First iteration of the average; generate variables
AvgIter1 = squeeze(mean(vidFrames,4));   % Double
Isum = 0*AvgIter1;
Nsum = 0*AvgIter1;
OneImage = Nsum +1;


% Second(final) iteration of the average
for i = 1:numFrames
    RawFrame = double(squeeze(vidFrames(:,:,1,i)));
    ValidPixels = double((RawFrame./AvgIter1 > LowCut) & (RawFrame./AvgIter1 < HiCut));
    Nsum = Nsum + ValidPixels;
    Isum = Isum + ValidPixels.*RawFrame;
end
AvgIter2 = Isum./Nsum;
fprintf('Average background calculated. '); toc


% Dynamic range readjustment; overwrites the loaded movie data to use less memory
for i = 1:numFrames
    RawFrame = double(squeeze(vidFrames(:,:,1,i)));
    ProcFrame = RawFrame./AvgIter2;                             % Normalize to the average
    ProcFrame = (1 - Exposures)*OneImage + Exposures*ProcFrame; % Linear intensity rescale for physical dark = 0
    ProcFrame(ProcFrame < 0) = 0;                               % Zero the negative values
    ProcFrame = ProcFrame*255*(1-HeadRoom*(1-LowCut));          % 8-bit range & prevent saturation
    vidFrames(:,:,1,i) = uint8(ProcFrame);                      % overwrite the frame
end
fprintf('Processing complete. Saving...'); toc


% Save the processed movie
v = VideoWriter(ProcName,'Grayscale AVI');
v.FrameRate = 15;
open(v)
for i = 1:numFrames
    frameM = vidFrames(:,:,1,i);
    writeVideo(v,frameM);
end
close(v)


% Generate the all-drop trace
CombinedDropPattern = AvgIter2 - AvgIter1;
CombinedDropPattern = max(max(CombinedDropPattern)) - CombinedDropPattern;

end
