% This script calculates fluorescence intensity and detects calcium events.
% Event detection is based on a clustering score derived from the brightest pixels in an expanded ROI.
% NOTE: Please modify file paths according to your local environment
% INPUT:
% - Motion-corrected calcium imaging data (filtered and unfiltered .tif)
% - Binary ROI masks (IC1.tif, IC2.tif, ...)
%
% OUTPUT:
% - Fluorescence intensity and event detection results (.csv)
% - ROI masks and related information (.mat, .tif)
%
% PARAMETERS:
% - nx: scaling factor for expanded ROI (default = 5)
% - percentage: fraction of brightest pixels used for clustering score (default = 20)
% - smooth: window size for temporal smoothing (default = 10)
% - thr: threshold for event detection (default = 0.4)
clearvars
mainfolder='C:\Users\Luchetti\Desktop\532p';
roifoldername='C:\Users\Luchetti\Desktop\532p\cell\common_ROI_eachday\intensity_res\ROIset\';
resfoldername = 'C:\Users\Luchetti\Desktop\532p\cell\common_ROI_eachday\intensity_res\';
mouse='#532';
year='2018';
month='11';
cd(mainfolder);
sdf='sd.csv';
sd_sheet=csvread(sdf);


for I=15:16
%v3: simplified, only 3 outputs in the xls:
% -frame number
% -smoothed score (default window 5)
% -spikes (default threshold 0.4)

%input:
%1. neuron acitvity video, named Input.tif
%2. IC roi binary, named IC1, IC2, ICX...
%3. results folder
%4. multiplication factor between small and big area
%5. percentage of brightest pixels for conf score
%6. sliding window for data smoothing of result
%7. index threshold for spike detection


clearvars -except I mainfolder roifoldername resfoldername month year sd_sheet mouse

num_day=sd_sheet(I,1);

if num_day<10
    day=['0',num2str(num_day)];
else
    day=num2str(num_day); 
end

switch sd_sheet(I,2)
    case 1
        session='a1';
    case 2
        session='a2';
    case 3
        session='b1';
    case 4
        session='b2';
    case 5
        session='c1';
    case 6
        session='c2';
end
date_name=[year,month,day];
FileName_input=[mainfolder,'\',date_name,'\',session];

%load the filtered video for clustering 
input_name = [year,month,day,mouse,session,'MCfilterDFFshifted.tif'];
fname=input_name;
%load the unfiltered video for intensity
input2_name =[year,month,day,mouse,session,'MCnonfilterDFFshifted.tif'];
fname2=input2_name;

%results folder
BaseName='session';
FileName1=[resfoldername,BaseName,num2str(I)];
mkdir(FileName1);
resfolder=FileName1;


BaseName='ROIs'; %changes folder to results folder
FileName2=[roifoldername,BaseName,num2str(I)];
roifolder = FileName2;
cd(roifolder);


%set how much larger the BIGROI area should be, default 5
nx = 5;

%choose percentage of pixels for the confidence score, default 20 
percentage = 20;

%sliding window for data smoothing of result
smooth = 10;

%index threshold for a spike to be detected
thr = 0.4;

%Note: how to prepare ROIs from the Mosaic output:
%1. open IC
%2. change to 8 bit
%3. set contrast bar to maximum
%4. make binary
%5. save
%(or use ImageJ macro: convert to 8bitcontra)

%=================%
%output:
%1. values of df local df and confidence score (.xls)
%2. BIGROI.mat
%3. SMALLROI.mat
%4. Picture of the ROI plus the two elipses (.jpeg)


cd(roifolder); %counts ROIs in the folder
roifiles = dir('IC*');
totroi= numel(roifiles);
frames = ceil(smooth/2);
cd(FileName_input);%change directory to input folder
info = imfinfo(fname);
num_images = numel(info);

totalroispikeCL=zeros(num_images,totroi);
totalroispikeFP=zeros(num_images,totroi);
totalroispikeDF=zeros(num_images,totroi);

for kroi = 1:totroi

%load the ROI
    X=[num2str(I),'_',month,day,session,'_IC',num2str(kroi)];
     display(X)

cd(roifolder);
BaseName='IC';
RoiName=[BaseName,num2str(kroi),'.tif']; %allow this to be visualized
bw = imread(RoiName);%load image of IC,regionfill?

%%identify ellipse
%use region props

SEprops = regionprops(bw, 'Orientation', 'MajorAxisLength', ...
    'MinorAxisLength', 'Eccentricity', 'Centroid'); % SE is small elipse

%Here are the measurements for the 255th object:
SEprops(255);

%representation/sanity check
imshow(bw)
hold on

phi = linspace(0,2*pi,50); % SE is small elipse
SEcosphi = cos(phi);
SEsinphi = sin(phi);


k = 255;
    SExbar = SEprops(k).Centroid(1);
    SEybar = SEprops(k).Centroid(2);
    
    SEa = SEprops(k).MajorAxisLength/2;
    SEb = SEprops(k).MinorAxisLength/2;
    
    SEtheta = pi*SEprops(k).Orientation/180;
    SER = [ cos(SEtheta)   sin(SEtheta)
         -sin(SEtheta)   cos(SEtheta)];
    
    SExy = [SEa*SEcosphi; SEb*SEsinphi];
    SExy = SER*SExy;
    
    SEx = SExy(1,:) + SExbar;
    SEy = SExy(2,:) + SEybar;

    plot(SEx,SEy,'r','LineWidth',2);

hold off
%turn those stats into an ellipse roi that is a mask
SMALLROI = roipoly(bw, SEx, SEy);

%%create a second bigger ellipse with nx times the Area
%ellipse area A = ? a b % nA ? sqrt(nx)a sqrt(nx)b 

phi = linspace(0,2*pi,50); % BE is big elipse
BEcosphi = cos(phi);
BEsinphi = sin(phi);


hold on

BExbar = SEprops(k).Centroid(1); 
    BEybar = SEprops(k).Centroid(2);
    
    BEa = SEprops(k).MajorAxisLength/2*sqrt(nx);
    BEb = SEprops(k).MinorAxisLength/2*sqrt(nx);
    
    BEtheta = pi*SEprops(k).Orientation/180;
    BER = [ cos(BEtheta)   sin(BEtheta)
         -sin(BEtheta)   cos(BEtheta)];
    
    BExy = [BEa*BEcosphi; BEb*BEsinphi];
    BExy = BER*BExy;
    
    BEx = BExy(1,:) + BExbar;
    BEy = BExy(2,:) + BEybar;

    plot(BEx,BEy,'r','LineWidth',2);

hold off

BIGROI = roipoly(bw, BEx, BEy); %turned into mask

%%save info of both ellipses as a structure


%based on http://blogs.mathworks.com/steve/2009/04/02/matlab-r2009a-imread-and-multipage-tiffs/
cd(FileName_input);%change directory to input folder
info = imfinfo(fname);
num_images = numel(info);

confmatrix = zeros (num_images,6);  %preallocate matrix for the final results
stats = regionprops (BIGROI,'pixellist'); %listing pixels of the two ROIs
x = stats.PixelList(:,1);
y = stats.PixelList(:,2);
statssmall = regionprops (SMALLROI,'pixellist');
xsmall = statssmall.PixelList(:,1);
ysmall = statssmall.PixelList(:,2);
avglist = zeros (num_images,4); % prepare matrix for average intensity for evry frame

for k = 1:num_images; % this FOR is the confidence value for every frame in the stack
    A = imread(fname, k, 'Info', info);
    roipixbig = impixel (A,x,y); %A is the image % found all pixels in the ROI and their intensity
    roipixsmall = impixel (A,xsmall,ysmall);
        %this calculates the DF between average intensity of small of big ROI
        avgbig = mean (roipixbig(:,1)); %average intensity bigroi
        avgsmall = mean (roipixsmall(:,1)); %average intensity smallroi
        avglist(k,1)=avgbig; % listing big ROI intensity for every frame
        avglist(k,2)=avgsmall;% listing big ROI intensity for every frame
        avglist(k,3)= avgsmall - avgbig;
        avg_int_allframe = mean(mean(A));
        avglist(k,4) = avgbig - avg_int_allframe;
        df = avgsmall - avgbig; %delta F
        dff = df / avgbig; %delta F / F
    pixeldata = [x y roipixbig(:,1)];
    sortedpixeldata = sortrows(pixeldata,-3); 
    top20a = size(sortedpixeldata); %finding the total amount of pixels in bigroi
    top20b = top20a(:,1); %finding the total amount of pixels in bigroi
    top20c = round(top20b/100*percentage);
    top20d = sortedpixeldata(1:top20c,:); % picked the X% brightest pixels
    x1 = top20d(:,1);
    y1 = top20d(:,2);
    mask = mat2gray (SMALLROI);
    findpixelsmall = impixel(mask,x1,y1); % multiplied top20 pixels by mask
    score = (sum(findpixelsmall(:,1)))/top20c; % final score is ratio of the X% brightest pixels from big area in small area
    confmatrix(k,1)=k; %first column is frame number
    confmatrix(k,2)=df; %second column is delta fluorescence
    confmatrix(k,3)=dff; %third column is delta fluorescence by local background fluorescence
    confmatrix(k,4)=score; %fourth column is conf value
end

    
%Calculating smoothed data for conf value

for k2 = 1:num_images;

a1 = (k2-frames);
     if a1 < 1 
     a1 = 1;
     end

a2 = (k2+frames);
     if a2 > num_images 
     a2 = num_images;
     end

smscore = mean(confmatrix(a1:a2,4));
confmatrix(k2,5)=smscore; %fifth column is smoothed conf value data    

end

%Detecting spikes by using the smoothed confindex and a threshold


for k3 = 1:num_images;
    
   if confmatrix(k3,5) > thr
   spiked = 1;
   else
   spiked = 0;
   end
   
confmatrix(k3,6)=spiked; %sixth column is spikes detected 
    
end


confmatrixv3(:,1)=confmatrix(:,1); %frames
confmatrixv3(:,2)=confmatrix(:,5); %smoother score
confmatrixv3(:,3)=confmatrix(:,6); %spikes

cd(FileName_input);%change directory to input folder
avglist2 = zeros (num_images,4); % prepare matrix for average intensity for unfiltered video
info2 = imfinfo(fname2);

for k = 1:num_images; % this FOR is the confidence value for every frame in the stack
    A2 = imread(fname2, k, 'Info', info2);
    roipixbig = impixel (A2,x,y); %A is the image % found all pixels in the ROI and their intensity
    roipixsmall = impixel (A2,xsmall,ysmall);
        %this calculates the DF between average intensity of small of big ROI
        avgbig = mean (roipixbig(:,1)); %average intensity bigroi
        avgsmall = mean (roipixsmall(:,1)); %average intensity smallroi
        avglist2(k,1)=avgbig; % listing big ROI intensity for every frame
        avglist2(k,2)=avgsmall;% listing big ROI intensity for every frame
        avglist2(k,3)= avgsmall - avgbig;
        avg_int_allframe = mean(mean(A2));
        avglist2(k,4) = avgsmall- avg_int_allframe;% listing small ROI-all pix intensity for every frame
        df = avgsmall - avgbig; %delta F
        dff = df / avgbig; %delta F / F
    
end

   fluo = avglist2(:,4); %small - all
    S = std(fluo);
SDdata = zeros(num_images,3);

    for k5 = 1:num_images;
        SDdata(k5,1)=avglist2(k5,4)./S; %(small - all)/SD
    end
        
  SD =  zeros(num_images,3);   
  SD(:,1) = avglist2(:,4); % substracted intensity small- all
  SD(:,2) = SDdata(:,1);
  
  for k6 = 1:num_images;
    if confmatrix(k6,5) > thr
   matchv = SD(k6,2);
   else
   matchv = 0;
    end
    
   SD(k6,3) = matchv;
   
  end
 
 SDspike = zeros(num_images,1); 
  for k7 = 1:num_images;
    if  SD(k7,3) > 1.54
   intspike = 1;
   else
   intspike = 0;
    end
    
  SDspike(k7,1) = intspike;
   
  end

  
finaldata = zeros(num_images,5);
finaldata(:,1) = avglist2(:,4); % substracted intensity small-all
finaldata(:,2) = SDdata(:,1);  %small-all/SD
finaldata(:,3) = confmatrixv3(:,2); % smoother score
finaldata(:,4) = SD(:,3);
finaldata(:,5) = SDspike(:,1);

totalroispikeCL(:,kroi)=SDspike(:,1);

%foopsi spike
[c,b,c1,g,sn,sp] = constrained_foopsi(SDdata(:,1));
foopsi=zeros(num_images,3);
foopsi(:,1)=SDdata(:,1);
foopsi(:,2)=sp;
for j=1:num_images
   if sp(j,1)>0.17
       foopsi(j,3)=1;
   end
end
totalroispikeFP(:,kroi)=foopsi(j,3);

%DF
DF=zeros(num_images,4);

Df(1,2)=0;
Y=diff(SDdata(:,1));
DF(:,1)=SDdata(:,1);
DF(2:num_images,2)=Y;
DF(:,3)=confmatrix(:,6);
for j=1:num_images
if j==1
    continue
end
   if Y(j-1,1)>0&&confmatrix(j,6)==1
    DF(j,4)=1;
   end

end
totalroispikeDF(:,kroi)=DF(j,4);

%=================%
%output:
%1. values of df, local df and confidence score (.xls)
%2. BIGROI.mat
%3. SMALLROI.mat
%4. Picture of the ROI plus the two elipses (.tif)


cd(resfolder);

BaseName='\IC_'; %changes folder to results folder
FileName=[resfolder,BaseName,num2str(kroi)];
mkdir(FileName);
resfoldersub = FileName;
cd(resfoldersub);

BaseName='IC_';
FileName=[BaseName,num2str(kroi),'.tif'];
saveas(gcf,FileName)

BaseName='bigroi_';
FileName=[BaseName,num2str(kroi)];
save(FileName, 'BIGROI');

BaseName='smallroi_';
FileName=[BaseName,num2str(kroi)];
save(FileName, 'SMALLROI');

BaseName='centroid_';
FileName=[BaseName,num2str(kroi)];
SEdata = SEprops(255);
save(FileName, 'SEdata');

% BaseName='results_';
% FileName=[BaseName,num2str(kroi),'.csv'];
% csvwrite(FileName,confmatrixv3);

BaseName='intensity_';
FileName=[BaseName,num2str(kroi),'.csv'];
csvwrite(FileName,finaldata);

BaseName='foopsi_';
FileName=[BaseName,num2str(kroi),'.csv'];
csvwrite(FileName,foopsi);

BaseName='DF_';
FileName=[BaseName,num2str(kroi),'.csv'];
csvwrite(FileName,DF);

% BaseName='SD_';
% FileName=[BaseName,num2str(kroi),'.csv'];
% csvwrite(FileName,SD);



cd ../ %moving out of results directory

end

filenameT='totalspikeCL.csv';
csvwrite(filenameT,totalroispikeCL);

filenameT='totalspikeFP.csv';
csvwrite(filenameT,totalroispikeFP);

filenameT='totalspikeDF.csv';
csvwrite(filenameT,totalroispikeDF);

end