preprocessH5.m

function preprocessH5(img_files,dark_files,ts_file,window_y,window_x,hot_pix_fcn,cycle_frame_num,save_folder,dark_save_folder,varargin)
%preprocessH5 Preprocess raw h5 files, doing dark offset and calculating
%windowed var and mean
%   Detailed explanation goes here
flag = false;
p = inputParser;
addOptional(p,'fileInd',1:numel(img_files));
addOptional(p,'temporal',false);
addOptional(p,'sectionFrameNum',1000);
parse(p,varargin{:});

fileInd = p.Results.fileInd;
temporal = p.Results.temporal;
sectionFrameNum = p.Results.sectionFrameNum;

frame_num_total = 0;
for k = fileInd
    frame_num_total = frame_num_total + ...
        numel(h5read(fullfile(img_files(k).folder,img_files(k).name),'/imgs',[1 1 1],[1 1 Inf]));
end

frame_num = 0;
frame_step = 100;

% timestamps
if isempty(ts_file)
    consider_skipped_frame = false;
else
    consider_skipped_frame = true;
end

if consider_skipped_frame
    ts = h5read(ts_file,'/timestamps');
    ts_diff = diff(ts);
    ts_diff_scale = double(ts_diff)./mean(double(ts_diff));
    ts_diff_scale = round(ts_diff_scale); % how many actual frames between each
end

save_file = 'dark_preprocessed.mat';
if ~exist(dark_save_folder,'dir')
    mkdir(dark_save_folder);
end

for dark_ind = 1:numel(dark_files)
    file_frame_num = h5info(fullfile(dark_files(dark_ind).folder,dark_files(dark_ind).name),'/imgs').Dataspace.Size(3);
    for section_ind = 1:floor(file_frame_num/frame_step)
        start_loc = [1 1 1+(section_ind-1)*frame_step];
        dark_imgs_file_section = h5read(fullfile(dark_files(dark_ind).folder,dark_files(dark_ind).name),...
            '/imgs',start_loc,[Inf Inf frame_step]);
        dark_imgs_file_section = permute(dark_imgs_file_section,[2 1 3]);

        frame_num = frame_num + frame_step;

        if section_ind == 1 && dark_ind == 1 % first file, first section
            dark_var_pix = var(double(dark_imgs_file_section),0,3)*frame_step;
            dark_mean_pix = mean(double(dark_imgs_file_section),3)*frame_step;
        else
            dark_var_pix = dark_var_pix + var(double(dark_imgs_file_section),0,3)*frame_step;
            dark_mean_pix = dark_mean_pix + mean(double(dark_imgs_file_section),3)*frame_step;
        end
    end
end
dark_var_pix = dark_var_pix./frame_num;
dark_mean_pix = dark_mean_pix./frame_num;

% get indices for each window
[Y,X] = meshgrid(1:size(dark_var_pix,1),1:size(dark_var_pix,2));
X = X'; Y = Y';
X = uint32(X); Y = uint32(Y);
x_coor_start = reshapeWindow(X,window_y,window_x); x_coor_start = x_coor_start(1,:);
y_coor_start = reshapeWindow(Y,window_y,window_x); y_coor_start = y_coor_start(1,:);

% window dark var and mean (pix x window)
dark_var_pix = reshapeWindow(dark_var_pix,window_y,window_x);
dark_mean_pix = reshapeWindow(dark_mean_pix,window_y,window_x);
window_num = size(dark_var_pix,2);

% find hot pixels
hot_pix = hot_pix_fcn(dark_mean_pix,dark_var_pix);

% calculate dark var window
dark_var_windowed = zeros(window_num,1);
dark_mean_windowed = zeros(window_num,1);
for window_ind = 1:window_num
    dark_var_windowed(window_ind) = mean(dark_var_pix(~hot_pix(:,window_ind),window_ind));
    dark_mean_windowed(window_ind) = mean(dark_mean_pix(~hot_pix(:,window_ind),window_ind));
end

% save
save(fullfile(dark_save_folder,save_file),'dark_var_windowed','dark_var_pix','dark_mean_windowed','dark_mean_pix','x_coor_start','y_coor_start');

if consider_skipped_frame
    % frame indices based on skipped time
    % first find what the true frame index is for each frame
    true_frame_ind = zeros(size(ts));
    true_frame_ind(1) = 1;
    for ts_diff_ind = 1:numel(ts_diff_scale)
        true_frame_ind(ts_diff_ind + 1) = true_frame_ind(ts_diff_ind) + ts_diff_scale(ts_diff_ind);
    end
    true_frame_ind = true_frame_ind(1:numel(ts));
    clear frame_cum

    % then map what frames should be saved to which saved frame ind
    save_frame_ind = zeros(size(ts));
    save_frame_ind(1) = 1;
    frame_cum = 1; % cumulative frame num

    % 1. figure out if there are skipped frames. If there are, shift the
    % save_frame_ind.
    % 2. figure out if the intensity time course matches what we expect
    % with the source detector separation. If not, shift the
    % save_frame_ind.

    % Step 1
    for ts_diff_ind = 1:numel(ts_diff_scale)
        if (ts_diff_scale(ts_diff_ind)) > cycle_frame_num && (ts_diff_ind == 1) % if there is large gap in ts at the beginning
            ts_diff_scale(ts_diff_ind) = 1; % assume thae first frame does not count
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum) + (1:ts_diff_scale(ts_diff_ind))';
            save_frame_ind = save_frame_ind + 1; 
            flag = true;
        elseif (ts_diff_scale(ts_diff_ind) > 1) && (ts_diff_ind > cycle_frame_num) % when there are skipped frames after the beginning
            % flag = true;
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum) + (1:ts_diff_scale(ts_diff_ind))';       % define the next frame normally (as if no skipped frames)
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)-1) = save_frame_ind((frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)-1) - cycle_frame_num); % since we are adding a filler frame without data, use the data from previous cycle.
            save_frame_ind(frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum+ts_diff_scale(ts_diff_ind)) - ts_diff_scale(ts_diff_ind) + 1;
        elseif (ts_diff_scale(ts_diff_ind) > 1) && (ts_diff_ind <= cycle_frame_num) % when there are skipped frames at the beginning
            % flag = true;
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum) + (1:ts_diff_scale(ts_diff_ind))';
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)-1) = save_frame_ind((frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)-1)) + cycle_frame_num - (ts_diff_scale(ts_diff_ind) - 1); % since we are adding a filler frame without data, use the data from previous cycle. This line adds a frame with the same illumination as the lost frame in the gap where the frame is lost. 
            save_frame_ind(frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum+ts_diff_scale(ts_diff_ind)) - ts_diff_scale(ts_diff_ind) + 1; % shifting back the frames such that e.g. og frame 2 is now used as saved frame 3
        else        
            save_frame_ind(frame_cum+1:frame_cum+ts_diff_scale(ts_diff_ind)) = save_frame_ind(frame_cum) + (1:ts_diff_scale(ts_diff_ind))';
        end
        frame_cum = frame_cum + ts_diff_scale(ts_diff_ind);
    end
    if flag % fixes the number of frames increasing
        save_frame_ind(numel(ts)) = save_frame_ind(numel(ts)-cycle_frame_num);
    end
    save_frame_ind = save_frame_ind(1:numel(ts));
else
    save_frame_ind = 1:frame_num_total;
    true_frame_ind = save_frame_ind;
end

% source ind of raw frames
source_ind_all = mod((1:numel(save_frame_ind)) - 1,cycle_frame_num);
source_ind_all = source_ind_all + 1;
raw_frame_source_ind = mod(true_frame_ind - 1,cycle_frame_num);
raw_frame_source_ind = raw_frame_source_ind + 1;

% sort files
file_names = {};
for f_ind = 1:numel(img_files)
    file_names{f_ind} = img_files(f_ind).name;
end
[~,I] = natsort(file_names);
img_files = img_files(I);

% illuminated files
file_frame_num = h5info(fullfile(img_files(1).folder,img_files(1).name),'/imgs').Dataspace.Size(3);
section_num = ceil(numel(save_frame_ind)./sectionFrameNum);

% Step 2
% define the expected intensity pattern based on sds

% get the actual intensity pattern

% % % load first section
% save_section_frame_ind = save_frame_ind(1:sectionFrameNum); % step 1 already applied
% imgs_section = getFrames(img_files,save_section_frame_ind,file_frame_num);
% imgs_section = permute(imgs_section,[2 1 3]);
% imgs_section = reshapeWindow(imgs_section,window_y,window_x);
% imgs_section_offset = double(imgs_section) - repmat(dark_mean_pix,[1 1 size(imgs_section,3)]);
% intensity_pattern = squeeze(mean(imgs_section_offset,[1 2]));
% intensity_pattern = mean(reshape(intensity_pattern(1:cycle_frame_num*10),cycle_frame_num,[]),2); % get the averaged pattern

% compare and identify any shifts


% shift the frames

if ~exist(save_folder,'dir')
    mkdir(save_folder);
end
for s_ind = 1:section_num
    disp(['Section # ' num2str(s_ind)]);
    save_file = ['section' num2str(s_ind) '_preprocessed.mat'];
    section_frame_ind = (s_ind-1)*sectionFrameNum + (1:sectionFrameNum);
    save_section_frame_ind = save_frame_ind(section_frame_ind);
    source_ind = source_ind_all(section_frame_ind);

    imgs_section = getFrames(img_files,save_section_frame_ind,file_frame_num);
    imgs_section = permute(imgs_section,[2 1 3]);
    imgs_section = reshapeWindow(imgs_section,window_y,window_x);

    % dark offset
    imgs_section_offset = double(imgs_section) - repmat(dark_mean_pix,[1 1 size(imgs_section,3)]);
    clear imgs_section

    % mean image
    for s = 1:cycle_frame_num
        img_mean(:,:,s) = mean(imgs_section_offset(:,:,source_ind == s),3);
        img_mean_frame_num(s) = sum(source_ind == s);
    end

    % window var and mean (pix x window x frame)
    var_windowed = zeros(window_num,size(imgs_section_offset,3));
    mean_windowed = zeros(window_num,size(imgs_section_offset,3));
    for window_ind = 1:window_num
        var_windowed(window_ind,:) = squeeze(var(imgs_section_offset(~hot_pix(:,window_ind),window_ind,:),0,1));
        mean_windowed(window_ind,:) = squeeze(mean(imgs_section_offset(~hot_pix(:,window_ind),window_ind,:),1));
    end
    if consider_skipped_frame
        save(fullfile(save_folder,save_file),...
            'source_ind','mean_windowed','var_windowed','img_mean','img_mean_frame_num','hot_pix','x_coor_start','y_coor_start','ts','save_frame_ind','raw_frame_source_ind','-v7.3');
    else
        save(fullfile(save_folder,save_file),...
            'source_ind','mean_windowed','var_windowed','img_mean','img_mean_frame_num','hot_pix','x_coor_start','y_coor_start','save_frame_ind','raw_frame_source_ind','-v7.3');
    end
end
end