%% Set up paths and make the initial file list
datadir = '/Users/mestapa/Documents/Research_Maine/NSF-STTR/MBARI_SES_data/';
filelist = dir(datadir);

%% Find the names of jpg files in the directory
rawfnams = [];
for i = 1:length(filelist)
    thisname = filelist(i).name;
    if length(thisname)>10 && strcmp(thisname(end-2:end),'jpg')
        rawfnams = [rawfnams;thisname];
    end
end

%% Read the image date-time from the filename
for i = 1:length(rawfnams)
    rawdates(i,:) = datevec(rawfnams(i,15:27),'yymmdd-HHMMSS');
end

%% Get rid of the second image in each pair (collected in the same hour)
% The first image is backlit, second is front-lit
testdates = datenum([rawdates(:,1:4),zeros(length(rawdates),2)]);
[imgdates,fileinds] = unique(testdates);
imgfnams = cellstr(rawfnams(fileinds,:));
clear raw* test* this* filelist i

%% Read the "red" channel of each image into a big matrix. CH added this note that 2219 and 2650 are pixel dimensions
allimg = zeros(2219,2650,length(imgfnams),'uint8');
for i = 1:length(imgfnams)
    if rem(i,100)==0;disp(num2str(i));end
    junk = imread(imgfnams{i});
    allimg(:,:,i) = junk(:,:,1);
    clear junk;
end
clear i;

%% Compute the variance in each pixel and find the image boundaries
for i = 1:size(allimg,1)
    thisslice = squeeze(single(allimg(i,:,:)));
    allimgvar(i,:) = var(thisslice,[],2)';
    if rem(i,100)==0;disp(num2str(i));end
end
edgethres = 100;
edgemask = ones(size(allimgvar));
edgemask(find(allimgvar<edgethres)) = NaN;
clear i thisslice;

%% Compute the median brightness in each image and find the brightest 10%
for i = 1:size(allimg,3)
    thisslice = single(allimg(:,:,i)).*edgemask;
    allmedbright(i,1) = median(thisslice(:),'omitnan');
    if rem(i,100)==0;disp(num2str(i));end
end

% prc95 = prctile(allmedbright,95);
% prc95ind = find(allmedbright>prc95);
% allgtprc95 = allimg(:,:,prc95ind);

prc90 = prctile(allmedbright,90);
prc90ind = find(allmedbright>prc90);
allgtprc90 = allimg(:,:,prc90ind);


%% Compute the background image as the mean of the brightest 10 pixel values from the 90th percentile images and then 
% smooth it with a median filter
rawbgr = mean(maxk(double(allgtprc90).*repmat(edgemask,[1 1 size(allgtprc90,3)]),10,3),3);
filtwin = [100 100];
bgr = medfilt2(rawbgr,filtwin);
clear i thisslice

%% Save the important output stuff
save SESdata.mat allimgvar allmedbright bgr edgemask edgethres filtwin imgdates imgfnams prc90 prc90ind rawbgr;