%% Script for analyzing 2p Ca-imaging series during imaging and stimulation

%Parameter setting
clear all; close all;
 
numROIs = 20;   % total number of ROIs;    
intrstm = 15;   % inter-stimulus time (s);
nSwps = 8;      % number of sweeps;
rf = 1346;      % id number of the reference img;
tsr = 1346;     % id number of the T-series;
option = 1;     % 1 for drawing, 0 for loading ROIs

%% Loading data
ref = imread([num2str(rf),'_reference.png']);     mask = imread([num2str(rf),'_mask.png']);
d = dir(['**/*TSeries*',num2str(tsr),'*.ome.tif*']);  frames = length(d);
xml = dir(['**/*TSeries*',num2str(tsr),'*.xml*']);    [fileID,~] = fopen(xml.name);
s = textscan(fileID,'%s','Delimiter','"');
fclose(fileID); frm_prd = find(~cellfun(@isempty,strfind(s{1},'framePeriod')));
dt = str2num(s{1}{frm_prd(1)+2}); 
for i = 1:frames
    tser(:,:,i) = imread(d(i).name);
end

%% Average projection and fluorescence timecourse, to locate stim epochs
proj = mean(tser,3);
fl = mean(mean(tser,2),1); fl = fl(:);  plot(fl);
[a,~] = ginput(2);  [m,frm] = max(fl(a(1):a(2)));
stim = round((a(1)+frm : intrstm/dt : a(1)+ frm + (nSwps-1)*intrstm/dt));

%% Show data and define ROIs for neurons and background
hFig1 = figure(1);   a=50;  b=550;
imshow(imfuse(ref(a:b,a:b),mask(a:b,a:b),'falsecolor','ColorChannels',[2 1 0]),[]);
set(hFig1, 'Position', [750 100 600 600])
hFig2 = figure(2);
imagesc(proj);
set(hFig2, 'Position', [10 50 700 700])

%Draw or load ROIs
for j = 1:numROIs
    if option == 1
        h = imfreehand(gca);    Roi = createMask(h);
        g = mean(getPosition(h)); text(g(1),g(2),['ROI',num2str(j)]);
        name = strcat('TSer_',d(1).name(45:48),'_ROI',num2str(j),'.mat');    save(name,'Roi'); 
    else
        name = strcat('TSer_',num2str(rf),'_ROI',num2str(j),'.mat');   Roi = load(name); Roi = Roi.Roi;
    end
    fluo(:,j) = ROIfluo(tser,Roi);
end

if option == 1
    saveas(gcf,['ROIs',num2str(rf),'.png']);
end

%% Analyze fluorescence traces
for j = 1:numROIs
   df_f(:,j) = deltaf_f(fluo(:,j));
end

win = floor(intrstm/dt)-30;
A = zeros(win,numROIs);
for i = 1:nSwps
    A = df_f(stim(i)-19:stim(i)+win-20,:) + A;
end
A = A./nSwps;   %%A(17:22,:) = NaN;

%plot results
hFig4 = figure(4);
set(hFig4, 'units','normalized','outerposition',[0.3 0 0.4 1]);
time = [0: dt: (size(A,1)-1)*dt];
for j = 1:numROIs
    subplot(2,1,1); plot(time, A(:,j));
    refline(0,3*std(A(1:15,j)));
    subplot(2,1,2); [am(j), ta(j)] = exp_fit(A(23:end,j),dt);
    sd(j) = std(A(end-20:end,j));
    w = waitforbuttonpress;
end

%save a table with results
Table=zeros(numROIs,3);
Table(:,1)=am';
Table(:,2)=sd';
Table(:,3)=ta';
dlmwrite('Table.txt',Table,'delimiter','\t','precision',4);
dlmwrite('df_f.txt',df_f(stim(1)-10:end,:),'delimiter','\t','precision',4);

