% This script performs multi-day cell tracking across sessions.
% Cells are matched based on centroid distance and spatial footprint similarity.
% Duplicate matches are resolved using a similarity-distance score.
% NOTE: Please modify file paths according to your local environment.


clearvars;

roifoldername = 'C:\Users\ab011\Documents\MATLAB\common_ROI_input\ROIs';
inputfolder = 'C:\Users\ab011\Documents\MATLAB\common_ROI_input';
matfolder = 'C:\Users\ab011\Documents\MATLAB';
resfolder = 'C:\Users\ab011\Documents\MATLAB\common_ROI_result';

%image size(pixel)
pix_x=270;
pix_y=360;

totalday=3;
totalsession=3; %total のsession数を入力
active_thr=1;%何セッション以上activeなものをとるか
dis_thr=6;%centroid distanceのthreshold
SI_thr=0.4;%SIのthreshold
 nse=1;

totroisheet=zeros(1,totalday);%roiの数を保存しておく
basename='XYtotalspike';
 cd(inputfolder);
for t=1:totalsession
    filename=[basename,num2str(t),'.csv'];
   
    AB= ['totalspike',num2str(t),'=csvread(filename);'];
   eval(AB);
    
end


for s=1:totalday %roiの数をそれぞれ数える
    roifolder=[roifoldername,num2str(s)];
cd(roifolder); %counts ROIs in the ROIfolder1
roifiles = dir('IC*');
totroisheet(1,s)= numel(roifiles);

end
maxroinum= max(totroisheet);
areasheet = zeros(maxroinum,1);

%読み込んだROI, active cell numberを保存する空シートを作る
for i2=1:totalday
   AB= ['ROIs_session',num2str(i2),'=zeros(pix_x,pix_y,totroisheet(1,i2));'];
   eval(AB);
   
   AB= ['active_session',num2str(i2),'=zeros(totroisheet(1,i2),1);'];
   eval(AB);
   
%    if i2==totalday
%       continue 
%    end
      AB=['dactive_cell_num',num2str(i2),'=zeros(totroisheet(1,i2),totalday-i2);'];
eval(AB);

   if i2==totalday
      AB=['dactive_cell_num',num2str(i2),'=zeros(totroisheet(1,i2),totalday-1);'];
     eval(AB);
   end

end

%ROIを読み込んで保存する
for i3=1:totalday
    totro=totroisheet(1,i3);
     roifolder=[roifoldername,num2str(i3)];
cd(roifolder);
    
    for i4=1:totro
    BaseName='IC';
RoiName=[BaseName,num2str(i4),'.tif']; 
bw = imread(RoiName);%load image of IC,regionfill?
AB=['ROIs_session',num2str(i3),'(:,:,i4)=bw;'];
    eval(AB);    
        
    end
    
end



for i4=1:totalday %一つ目のセッション

     
    for k=1:totroisheet(1,i4)
         X=['Calculating distance and SI in session',num2str(i4),'cell#',num2str(k)];
     display(X)
       %ROI ひとつめ
AB=['bw1=ROIs_session',num2str(i4),'(:,:,k);'];
    eval(AB);    

Centroid = regionprops(bw1, 'Centroid'); 
c = struct2cell(Centroid);
cen1 =c{1};
          

for i5=1:totalday %二つ目のセッション
    if i4==totalday
        if i5==i4
            continue
        end
    else
        if i5<=i4
            continue
        end
    end

    check = zeros(totroisheet(1,i5),1);
    this_samecell=0;  
    
for kroi = 1:totroisheet(1,i5)
% if check(kroi,1)==1
%    continue 
% end

%load the ROI
AB=['bw2=ROIs_session',num2str(i5),'(:,:,kroi);'];
    eval(AB);    


Centroid = regionprops(bw2, 'Centroid'); 
c = struct2cell(Centroid);
cen2 =c{1};
a=regionprops(bw2, 'Area'); 
        ar = struct2cell(a);
        area2 =ar{255};

 D = pdist2(cen1,cen2);
 
 % SI=dot(reshape(bw1,pix_y*pix_x,1),reshape(bw2,pix_y*pix_x,1)/norm(reshape(bw1,pix_y*pix_x,1))/norm(reshape(bw2,pix_y*pix_x,1)));
  v1 = reshape(bw1, pix_y*pix_x, 1);
v2 = reshape(bw2, pix_y*pix_x, 1);
SI = dot(v1, v2) / (norm(v1) * norm(v2));

   if SI>SI_thr&&D<dis_thr %bw1とbw2が同じと判断された場合bw2のROI番号をactive_cell_numに登録
        this_samecell=cat(1,this_samecell,kroi); 

        
   end


end%kroi のend
        this_samecell(1,:)=[];

        
        if i4==totalday
            AB=['active_cell_num',num2str(i4),'{k,i5}=this_samecell;'];
            eval(AB);
        else
            AB=['active_cell_num',num2str(i4),'{k,i5-i4}=this_samecell;'];
            eval(AB);
            
        end
        
end%ふたつめのsession　end

    end
end

%%
for AC=1:totalday
    eval(['thisdata=active_cell_num',num2str(AC),';']);
      [R T]=size(thisdata);
    for AD=1:T%レーン
        thislane=thisdata(:,AD);
        if AC==totalday
           laneses=AD;
        else
             laneses=AD+AC;
        end
        for k=1:totroisheet(laneses)
             check_lane=zeros(length(thislane),1);
            for l=1:length(thislane)
                data=thislane{l};
                finds=find(data==k);
                if isempty(finds)==0
                    check_lane(l)=length(data);
                end
                

            end
            ch=sum(check_lane>0);
            if ch>1
                finish=0;
                for m=1:max(check_lane)
                    findm=find(check_lane==m);
                    if findm==0
                        continue
                    end
                    if finish==1
                        continue
                    end
               
                        AB=['bw1=ROIs_session',num2str(laneses),'(:,:,k);'];
                        eval(AB);
                        
                        Centroid = regionprops(bw2, 'Centroid');
                        c = struct2cell(Centroid);
                        cen1 =c{1};
                        a=regionprops(bw2, 'Area');
                        ar = struct2cell(a);
                        area2 =ar{1};
                        
                        sid=zeros(length(findm),1);
                        for si=1:length(findm)
                            AB=['bw2=ROIs_session',num2str(AC),'(:,:,findm(si));'];
                            eval(AB);
                            
                            Centroid = regionprops(bw2, 'Centroid');
                            c = struct2cell(Centroid);
                            cen2 =c{1};
                            a=regionprops(bw2, 'Area');
                            ar = struct2cell(a);
                            area2 =ar{1}; 
                            D = pdist2(cen1,cen2);
                            SI=dot(reshape(bw1,pix_y*pix_x,1),reshape(bw2,pix_y*pix_x,1)/norm(reshape(bw1,pix_y*pix_x,1))/norm(reshape(bw2,pix_y*pix_x,1)));
                            v1 = reshape(bw1, pix_y*pix_x, 1);
                            v2 = reshape(bw2, pix_y*pix_x, 1);
                            SI = dot(v1, v2) / (norm(v1) * norm(v2));
                            if D == 0
                                sid(si) = SI;
                            else
                                sid(si) = SI / D;
                            end
                            
                        end
                       maxsid=find(sid==max(sid));
                       champion=findm(maxsid);
                       
                       %championいがい、けす
                       for l2=1:length(check_lane)
                           if l2==champion
                               continue
                           end
                           if check_lane(l2)>0
                               thecande=thislane{l2};
                               thecande(thecande==k)=[];
                               if isempty(thecande)==1
                                   thecande=0;
                               end
                               thisdata{l2,AD}= thecande;
                           end
                           
                       end %end of l2
                       
                       finish=1;
                       
                end %end of m
 
            end %if ch
        end

    end
   eval(['active_cell_num',num2str(AC),'=thisdata;']); 
    
end

%%

%ダブリを消す
for e=1:totalday-1
         X=['DEL1_',num2str(e)];
     display(X)
     
    AB=['BB=active_cell_num',num2str(e),';'];
    eval(AB);
    
    for e2=1:totalday-1%消しに行くセッション
        if e2<=e
            continue
        end
        
        
        I=e2-e;%e2はBBのどの列
        
        for e3=1:totroisheet(1,e)
            thecan=BB{e3,I};
            if isempty(thecan)==1
                continue
            end
            
            if length(thecan)>1
                res=zeros(length(thecan),1);
                for kr=1:length(thecan)
                    num_c=thecan(kr);
                    AB=['thisdata=active_cell_num',num2str(e2),';'];
                     eval(AB);
                    [R T]=size(thisdata);
                    for t=1:T
                        t_ses=e2+t;
                        nc=thisdata{num_c,t};
                        if isempty(nc)==1
                            continue
                        end
                         target=t_ses-e;
                         num=BB{e3,target};
                        if isempty(num)==1
                            continue
                        end
                         if isempty(find(nc==num(1,1))==0)
                             res(kr)=res(kr)+1;
                         end
                    end
                    
                end %end of kr
                champcell= thecan(find(res==max(res)));
                ab=thecan==champcell(1,1);
                thecan(ab==0)=[];

            end
            BB{e3,I}=thecan;        
        end%end of e3
        
    end
            AB=['active_cell_num',num2str(e),'=BB;'];
        eval(AB);
end

%%
for e=1:totalday-1
    
     X=['DEL2_',num2str(e)];
     display(X)
    AB=['BB=active_cell_num',num2str(e),';'];
    eval(AB);
    
    e2=totalday;
    I=e2-e;
    
    for e3=1:totroisheet(1,e)
        thecan=BB{e3,I};
        if isempty(thecan)==1
            continue
        end
        
        if length(thecan)>1
            res=zeros(length(thecan),1);
            for kr=1:length(thecan)
                num_c=thecan(kr);
                AB=['thisdata=active_cell_num',num2str(e2),';'];
                eval(AB);
               
                    t=e;
                    nc=thisdata{num_c,t};
                    if isempty(nc)==1
                        continue
                    end

                    if isempty(find(e3==nc))==0
                        res(kr)=res(kr)+1;
                    end
       
                
            end %end of kr
            champcell= thecan(find(res==max(res)));
            ab=thecan==champcell;
            thecan(ab==0)=[];
             BB{e3,I}=thecan;
        end
      
             
    end
    
            AB=['active_cell_num',num2str(e),'=BB;'];
        eval(AB);
    
end

%%
% for AE=1:totalday-1
%     AB=['BB=active_cell_num',num2str(AE),';'];
%     eval(AB);
%     AB=['thisdata=active_cell_num',num2str(totalday),';'];
%     eval(AB);
%      I=totalday-AE;
%     for kr=1:totroisheet(1,totalday)
%        num=thisdata{kr,AE};
%        if num==0
%               continue
%        end 
%        if isempty(num)==1
%               continue
%        end
%        target=BB{num,I};
%        
%        if isempty(target)==1
%            
%            BB{num,I}=kr;
%        end
%     end 
%     AB=['active_cell_num',num2str(AE),'=BB;'];
%     eval(AB);
% 
% end

%%

for i4=1:totalday-1 %一つ目のセッション
             X=['DEL3_',num2str(e)];
     display(X)
    AB=['BB=active_cell_num',num2str(i4),';'];
    eval(AB);
    
    [R T]=size(BB);
    
    for i5=1:T
        t_ses=i4+i5;
        for k=1:R
            num=BB{k,i5};
            if length(num)==1
                AB=['dactive_cell_num',num2str(i4),'(k,i5)=num;'];
                eval(AB);
                
            end
            
            if length(num)>1
                %ROI ひとつめ
                AB=['bw1=ROIs_session',num2str(i4),'(:,:,k);'];
                eval(AB);
                
                Centroid = regionprops(bw1, 'Centroid');
                c = struct2cell(Centroid);
                cen1 =c{1};
                
                res=zeros(length(num),1);
                for kr=1:length(num)
                    this_cell=num(kr);
                    AB=['bw2=ROIs_session',num2str(t_ses),'(:,:, this_cell);'];
                    eval(AB);
                    Centroid = regionprops(bw2, 'Centroid');
                    c = struct2cell(Centroid);
                    cen2 =c{1};
                    
                    D = pdist2(cen1,cen2);
                   SI=dot(reshape(bw1,pix_y*pix_x,1),reshape(bw2,pix_y*pix_x,1)/norm(reshape(bw1,pix_y*pix_x,1))/norm(reshape(bw2,pix_y*pix_x,1)));
                    res(kr)=SI/D;
                end
                 champcell= num(find(res==max(res)));
                 ab=num==champcell;
                 num(ab==0)=[];
                  BB{k,i5}=num;
                  AB=['dactive_cell_num',num2str(i4),'(k,i5)=num;'];
                  eval(AB);
            end
        end%end of k
    end
     AB=['active_cell_num',num2str(i4),'=BB;'];
    eval(AB);
end
 

%%

%ダブリを消す
for e=1:totalday-1
         X=['DEL4_',num2str(e)];
     display(X)
     
    AB=['BB=dactive_cell_num',num2str(e),';'];
    eval(AB);
    
    for e2=1:totalday%消しに行くセッション
        if e2<=e
            continue
        end
        
        
        I=e2-e;%e2はBBのどの列
        
        for e3=1:totroisheet(1,e)
        num_c=BB(e3,I);
        if num_c==0
            continue
        end
        
      
        AB=['dactive_cell_num',num2str(e2),'(num_c,:)=0;'];
        eval(AB);

        end
        
    end
    
end
%%
for e=1:totalday-1
         X=['DEL5_',num2str(e)];
     display(X)
     
    AB=['BB=dactive_cell_num',num2str(e),';'];
    eval(AB);
    
    for e2=1:totalday%消しに行くセッション
        if e2<=e
            continue
        end
        
        
        I=e2-e;%e2はBBのどの列
        
        for e3=1:totroisheet(1,e)
        num_c=BB(e3,I);
        if num_c==0
            continue
        end
        for y=e+1:e2-1
            AB=['thisdata=dactive_cell_num',num2str(y),';'];
            eval(AB);
             I2=e2-y;%e2はどの列
             a=find(thisdata(:,I2)==num_c);
             if isempty(a)==0
                 for x=1:length(a)
                     I3=y-e;
                     if BB(e3,I3)==0
                         BB(e3,I3)=a(x);
                     end
                 end
                 thisdata(a(x),:)=0;
             end
             AB=['dactive_cell_num',num2str(y),'=thisdata;'];
             eval(AB);
        end
      
        end
        
    end
         AB=['dactive_cell_num',num2str(e),'=BB;'];
        eval(AB);
end


for AC=1:totalday
      X=['DEL6_',num2str(AC)];
     display(X)
     
    eval(['thisdata=dactive_cell_num',num2str(AC),';']);
      [R T]=size(thisdata);
    for AD=1:T%レーン
        thislane=thisdata(:,AD);
        if AC==totalday
           laneses=AD;
        else
             laneses=AD+AC;
        end
        for k=1:totroisheet(laneses)
             check_lane=zeros(length(thislane),1);
            for l=1:length(thislane)
                data=thislane(l);
                finds=find(data==k);
                if isempty(finds)==0
                    check_lane(l)=length(data);
                end
                

            end
            ch=sum(check_lane>0);
            if ch>1
                finish=0;
                for m=1:max(check_lane)
                    findm=find(check_lane==m);
                    if findm==0
                        continue
                    end
                    if finish==1
                        continue
                    end
               
                        AB=['bw1=ROIs_session',num2str(laneses),'(:,:,k);'];
                        eval(AB);
                        
                        Centroid = regionprops(bw2, 'Centroid');
                        c = struct2cell(Centroid);
                        cen1 =c{1};
                        a=regionprops(bw2, 'Area');
                        ar = struct2cell(a);
                        area2 =ar{1};
                        
                        sid=zeros(length(findm),1);
                        for si=1:length(findm)
                            AB=['bw2=ROIs_session',num2str(AC),'(:,:,findm(si));'];
                            eval(AB);
                            
                            Centroid = regionprops(bw2, 'Centroid');
                            c = struct2cell(Centroid);
                            cen2 =c{1};
                            a=regionprops(bw2, 'Area');
                            ar = struct2cell(a);
                            area2 =ar{1}; 
                            D = pdist2(cen1,cen2);
                            SI=dot(reshape(bw1,pix_y*pix_x,1),reshape(bw2,pix_y*pix_x,1)/norm(reshape(bw1,pix_y*pix_x,1))/norm(reshape(bw2,pix_y*pix_x,1)));
                        sid(si)=SI/D;
                        
                        end
                       maxsid=find(sid==max(sid));
                       champion=findm(maxsid);
                       
                       %championいがい、けす
                       for l2=1:length(check_lane)
                           if l2==champion
                               continue
                           end
                           if check_lane(l2)>0
                               thecande=thislane(l2);
                               thecande(thecande==k)=[];
                               if isempty(thecande)==1
                                   thecande=0;
                               end
                               thisdata(l2,AD)= thecande;
                           end
                           
                       end %end of l2
                       
                       finish=1;
                       
                end %end of m
 
            end %if ch
        end

    end
   eval(['dactive_cell_num',num2str(AC),'=thisdata;']); 
    
end


%%
%active_session
for t=1:totalday-1
     X=['Calculating active number in session',num2str(t)];
     display(X)
    
    AB=['active_cell_num=dactive_cell_num',num2str(t),';'];
    eval(AB);
    [m,n]=size(active_cell_num);
    active_cell_num(active_cell_num>0)=1;
    
    AB=['active_session',num2str(t),'=sum(active_cell_num,2);'];
    eval(AB);
end
%%
 
CC=0;
for j=1:totalday
     AB=['act=active_session',num2str(j),';'];
    eval(AB);
    
   for j2=1:totroisheet(1,j)
       
       if act(j2,1)>=active_thr
       CC=CC+1;
       ori_info(CC,1)=j;%original session
       ori_info(CC,2)=j2;%original cell number
       AB=['roi=ROIs_session',num2str(j),'(:,:,j2);'];
       eval(AB);
       
       total_ROIs(:,:,CC)=roi;
           
       end
          
   end
    
end


[a b]=size(ori_info);
final_ROInum=a;
del_check=zeros(final_ROInum,1);
for k=1:final_ROInum
    if  del_check(k)==1
        continue
    end
    this_ses=ori_info(k,1);
    this_roi=ori_info(k,2);
  
    for k2=1:final_ROInum
        if k==k2
            continue
        end
        
        if this_ses==ori_info(k2,1)
            if this_roi==ori_info(k2,2)
                del_check(k2)=1;
            end
        end
    end
end
 ori_info(del_check==1,:)=[];
[a b]=size(ori_info);
final_ROInum=a;
cell_numsheet=zeros(final_ROInum,totalday);

for m=1:totalday
        
       for m2=1:final_ROInum
           if ori_info(m2,1)==m
              cell_numsheet(m2,m)=ori_info(m2,2);
              
           else
               snum=ori_info(m2,1);
              if snum>m
              continue
              else
               AB=['data=dactive_cell_num',num2str(snum),';'];
               eval(AB);
               cell_numsheet(m2,m)=data(ori_info(m2,2),m-snum);
              end  
           end 
       end
end
   

 

%totalspike書き出し
cd(resfolder);
session=0;
for day=1:totalday

%     if day>9
%         nse=1;
%     else
%         nse=1;
%     end
  
for i=1:nse
    session=session+1;
   AB=['spikedata=totalspike',num2str(session),';'];
   eval(AB);
   A=spikedata(:,1:3);
   spikedata(:,1:3) = [];
 [a b]=size(spikedata);  
   nframe=a;
    csvdata=zeros(nframe,final_ROInum);
    cd(resfolder);
    for kroi=1:final_ROInum
      if cell_numsheet(kroi,day)>0
          csvdata(:,kroi)=spikedata(:,cell_numsheet(kroi,day));
          
      end
    end
    C = horzcat(A,csvdata);
      fileNaMe=['XYtotalspike',num2str(session),'.csv'];
  csvwrite(fileNaMe,C);
end

end

filename='workspace.mat';
save(filename);

imagesc(sum(total_ROIs,3))
mname ='merged.tif';
saveas(figure(1),mname);

fileNaMes='cell_numsheet.csv';
  csvwrite(fileNaMes,cell_numsheet);