% This script detects candidate cell ROIs from calcium imaging data.
% ROIs are extracted by frame-wise thresholding, size/shape filtering,
% and duplicate removal based on spatial similarity.
% NOTE: Please modify file paths according to your local environment.
%
% INPUT:
% - Frame-by-frame calcium imaging .tif images
% - Binary mask image (mask.tif)
% - Session information file (sd.csv)
%
% OUTPUT:
% - Binary ROI masks (IC1.tif, IC2.tif, ...)
% - Merged ROI image (merged.tif)
% - MATLAB workspace file (workspace.mat)
%
% PARAMETERS:
% - SD: threshold multiplier for image binarization
% - s: allowed ROI size range in pixels
% - spatial similarity threshold: used for duplicate ROI removal

clearvars 
mainfolder='F:\cell_detection\173ex';
resfolder = 'F:\cell_detection\173ex\result';
cd(mainfolder);
sdf='sd.csv';
sd_sheet=csvread(sdf);

I = imread('mask.tif');
I2 = im2double(I);
mask = I2;


for x=11:18
    
%     if j<7
%         cd(mainfolder);
%         I = imread('mask.tif');
%         I2 = im2double(I);
%         mask = I2;
%     else
%         cd(mainfolder);
%         I = imread('mask2.tif');
%         I2 = im2double(I);
%         mask = I2;
%         
%     end
    
clearvars -except x sd_sheet mainfolder  month resfolder mask
close all

matfolder='C:\Users\ab011\Documents\MATLAB';

x_px=320;
y_px=200;
num_month=sd_sheet(x,1);
if num_month<10
    month=['0',num2str(num_month)];
else
    month=num2str(num_month); 
end

num_day=sd_sheet(x,2);

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

switch sd_sheet(x,3)
    case 1
        session='a1';
    case 2
        session='a2';
    case 3
        session='b';
    case 4
        session='c1';
    case 5
        session='c2';
    case 6
        session='c2';
end

%Note: images should be 8 bit.
%first, make images of each frame
imagefolder=[mainfolder,'\',num2str(month),day,'\',session];
SD=sd_sheet(x,4);
if SD==0
    continue
    
end

s=[30 180]; %size of cell default (50 150)
cd(imagefolder);
imlist=dir('*.tif');
ggg=1000;
filter2=[];
filter2(:,:,1)=ones(y_px,x_px);

num_rois=0;
im_data=[];
gixel_size=[];

%create a mask
   display(strcat(month,day,session,'_SD=',num2str(SD)))
cd(imagefolder);


for i=1:numel(imlist); %number of images
    
    im=im2double(imread(imlist(i).name));
    im_data(:,:,i)=im;
    im_mean=mean(reshape(im,size(im,1)*size(im,2),1));
    im=im-im_mean; %全体の明るさから平均を?ｷし引く
    im_thr=mean(reshape(im,size(im,1)*size(im,2),1))+SD*std(reshape(im,size(im,1)*size(im,2),1));%3*std
    p=find(im>im_thr); 
    k=find(im<=im_thr);
    im(p)=1;
    im(k)=0;
    
   im=im2double(im); 
   im = im.*mask; %mask
    
    im = imfill(im, 'holes');
    labeledImage = bwlabel(im, 8);     
    coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
    %imagesc(coloredLabels);
    blobMeasurements = regionprops(labeledImage, im, 'all');   
    numberOfBlobs = size(blobMeasurements, 1);
    pixels=zeros(size(im,1),size(im,2), numberOfBlobs);
    ratx=[];
    raty=[];
    for k = 1 : numberOfBlobs           
        blobArea= blobMeasurements(k).Area;
        blobCentroid= blobMeasurements(k).Centroid;
        blobRoi= blobMeasurements(k).PixelList;
        
        pix=size(blobRoi);
        xx=blobRoi(:,2);
        xw=max(xx)-min(xx);
        
        yy=blobRoi(:,1);
        yw=max(yy)-min(yy);
        xy=xw/yw;
        yx=yw/xw;
        for jq=1:pix(1)
            pixels(xx(jq),yy(jq),k)=1;
        end
        pixels(:,:,k)= imfill(pixels(:,:,k), 'holes');
        ratx(k)=xy;
        raty(k)=yx;
     
    end
    
    allBlobAreas = [blobMeasurements.Area];
    allblobCentroid= [blobMeasurements.Centroid];
    allowableAreaIndexes = (allBlobAreas > s(1))&(allBlobAreas < s(2)& (ratx<2)&(raty<2));
    %allowableAreaIndexes=intersect(allowableAreaIndexes1,allowableAreaIndexes2);
    keeperIndexes = find(allowableAreaIndexes);
    
    keeperBlobsImage = ismember(labeledImage, keeperIndexes);
    gixels=pixels(:,:,keeperIndexes);
    gixel_size(i)=size(gixels,3);
    if isempty(gixels)
        num_rois=num_rois;
    else
        num_rois=num_rois+size(gixels,3);
        filter2=cat(3,filter2,gixels);
    end
   display(strcat(num2str(num_rois),' rois in ',num2str(i),' images'))
   
    
%     imagesc(sum(filter2,3))
%     drawnow
%     
    
end
sig=[];
filter2(:,:,1)=[];

 pass_first=ones(size(filter2,3),size(filter2,3));
pass_first=pass_first(:,1);

for i=1:size(filter2,3)
    
    display(strcat('first deletion',num2str(sum(pass_first))))
    if pass_first(i)==1
        for j=1:size(filter2,3)
            if i==j||pass_first(j)==0
            else
                 SI=dot(reshape(filter2(:,:,i),y_px*x_px,1),reshape(filter2(:,:,j),y_px*x_px,1))/norm(reshape(filter2(:,:,i),y_px*x_px,1))/norm(reshape(filter2(:,:,j),y_px*x_px,1));
%                 corr=corrcoef(sig(i,:),sig(j,:));
%                 corr=corr(2,1);
                if (SI>0.2) %Tempral相関の閾値0.9
                    if sum(sum(filter2(:,:,i)))>=sum(sum(filter2(:,:,j)))
                        pass_first(i)=0;
                    else
                        pass_first(j)=0;
                    end
                end
            end
        end
    else
    end
end


filter2(:,:,pass_first==0)=[];
imagesc(sum(filter2,3))


cd(resfolder);
BaseName='\output_'; %changes folder to results folder
FileName=[resfolder,BaseName,num2str(month),num2str(day),session,'_',num2str(SD),'sd'];
mkdir(FileName);
resfoldersub=FileName;
cd(resfoldersub);
filename='workspace.mat';
save(filename);
num_finalroi = length(filter2(1,1,:));

mname ='merged.tif';
saveas(figure(1),mname);

for i=1:num_finalroi; %save each ROIs as TIF
    
fig_data = filter2(:,:,i);
BaseName='IC';
FileName=[BaseName,num2str(i),'.tif'];
imwrite(fig_data,FileName);

end


end

cd(matfolder);