msrd_spatial_gitHub.m

%% Compute MSRD values for individual cells
clearvars; 

%% Specify directory path and file name
idr='/Users/input_directory/';

sub_dir=dir(sprintf('%s*data_type*',idr));
cs_all=cell(size(sub_dir,1),1);
for ii=1:size(sub_dir,1)
    cs_all{ii}=sub_dir(ii).name;
end

clear sub_dir;

for csc=1:numel(cs_all)
    csNm=cs_all{csc};
    
    %% Run code
    % Import position data
    pDtTb=readtable(sprintf('%s%s/%s_fulltb_position.csv',idr,...
        csNm,csNm),'HeaderLines',3);
    pDt=zeros(size(pDtTb,1),6);
    pDt(:,1:5)=[table2array(pDtTb(:,1:3)),table2array(pDtTb(:,7)),...
        table2array(pDtTb(:,8))-10^9+1];
    
    % Import position reference points
    pRfTb=readtable(sprintf('%s%s/%s_fulltb_reference_point.csv',...
        idr,csNm,csNm),'HeaderLines',3);
    pRf=table2array(pRfTb(:,1:3));
    
    % Import dorsal-ventral boundary
    pdvbTb=readtable(sprintf('%s%s/%s_fulltb_dv_boundary.csv',...
        idr,csNm,csNm),'HeaderLines',3);
    pdvb=table2array(pdvbTb(:,1:3));
    
    % Import epithelial-dorsal boundary
    pedbTb=readtable(sprintf('%s%s/%s_fulltb_ed_boundary.csv',...
        idr,csNm,csNm),'HeaderLines',3);
    pedb=table2array(pedbTb(:,1:3));
    
    clear pDtTb pRfTb pdvbTb pedbTb;
    
    % Specify origin and length scale
    origin=[pRf(2,1),pRf(2,2),min(pDt(:,3))];
    len_scale=110;
    
    % Interpolation function of boundary surfaces
    dvbSrf=scatteredInterpolant(pdvb(:,1),pdvb(:,2),pdvb(:,3));
    edbSrf=scatteredInterpolant(pedb(:,1),pedb(:,2),pedb(:,3));
    
    dv_cri=dvbSrf(pDt(:,1),pDt(:,2))-pDt(:,3);
    ed_cri=edbSrf(pDt(:,1),pDt(:,2))-pDt(:,3);
    
    % delete all epithelial cell data
    vn_id=unique(pDt(dv_cri<0,5));
    ds_id=unique(pDt((dv_cri>=0 & ed_cri<0),5));
    ep_id=unique(pDt(ed_cri>=0,5));
    ep_sol_id=setdiff(ep_id,ds_id);
    ep_sol_id=setdiff(ep_sol_id,vn_id);
    
    for ii=1:size(ep_sol_id)
        pDt=pDt(pDt(:,5)~=ep_sol_id(ii),:);
    end
    
    dv_cri=dvbSrf(pDt(:,1),pDt(:,2))-pDt(:,3);
    
    % specify region number for each data. 1:dorsal, 2:dorsal MPZ, 3:ventral
    % MPZ, 4:PSM.
    pDt(dv_cri<0,6)=4;
    pDt(dv_cri>=0,6)=1;
    pDt((pDt(:,2)>=pRf(1,2) & dv_cri<0),6)=3;
    pDt((pDt(:,2)>=pRf(1,2) & dv_cri>=0),6)=2;
    
    % Total cell number and max time
    vId=unique(pDt(:,5));
    ttNm=size(vId,1);
    
    % Create sorted array in terms of particle ID
    pPc=cell(ttNm,1);
    
    for ii=1:ttNm
        pPc{ii}=pDt(pDt(:,5)==vId(ii),:);
    end
    
    % Check all cells that crosses between region1 and region4
    chk=zeros(ttNm,1);
    for ii=1:ttNm
        reg_dff=abs(pPc{ii}(1:end-1,6)-pPc{ii}(2:end,6));
        reg_dff_cs=unique(reg_dff);
        if max(reg_dff)==3
            chk(ii)=1;
        end
    end
    
    % New data set that excludes crossing cells
    pPc=pPc(chk==0);
    vId=vId(chk==0);
    ttNm=size(vId,1);
    
    pDt=cell2mat(pPc);
    
    % Created sorted array in terms of time
    ttTm=max(pDt(:,4));
    pTm=cell(ttTm,1);
    for ii=1:ttTm
        pTm{ii}=pDt(pDt(:,4)==ii,:);
    end
    % 
    % % Compute fluidity index for each cell
    msrdDt=zeros(ttNm*ttTm,7);
    
    [dsc,neiMx,tmScl]=deal(1,8,15);
    
    tic;
    for ii=1:ttNm
        for jj=1:size(pPc{ii},1)-tmScl
            cenDt_i=pPc{ii}(jj,:);
            tm_i=cenDt_i(4);
            pps=find(pTm{tm_i}(:,5)==vId(ii));
            ds=sqrt((pTm{tm_i}(:,1)-pTm{tm_i}(pps,1)).^2+...
                (pTm{tm_i}(:,2)-pTm{tm_i}(pps,2)).^2+...
                (pTm{tm_i}(:,3)-pTm{tm_i}(pps,3)).^2);
            [~,indx]=sort(ds);
            nei_idx=zeros(neiMx,1);
            nei_chk=1;
            nei_cnt=2;
            while nei_chk<=8
                nei_id=indx(nei_cnt);
                if abs(pTm{tm_i}(pps,6)-pTm{tm_i}(nei_id,6))<2
                    nei_idx(nei_chk)=nei_id;
                    nei_chk=nei_chk+1;
                end            
                nei_cnt=nei_cnt+1;
            end
            neiId=pTm{tm_i}(nei_idx,5);            
    
            cenDt_f=pPc{ii}(jj+tmScl,:);
            tm_f=cenDt_f(4);
            msrdTmp=zeros(neiMx,1);
            msc=1;
            for kk=1:neiMx
                neiDt=pPc{vId==neiId(kk)};
                neiDt_i=neiDt(neiDt(:,4)==tm_i,:);
                neiDt_f=neiDt(neiDt(:,4)==tm_f,:);
                if isempty(neiDt_f)==0
                    msrdTmp(msc)=(((cenDt_f(1)-neiDt_f(1))...
                        -(cenDt_i(1)-neiDt_i(1)))^2+...
                        ((cenDt_f(2)-neiDt_f(2))...
                        -(cenDt_i(2)-neiDt_i(2)))^2+...
                        ((cenDt_f(3)-neiDt_f(3))...
                        -(cenDt_i(3)-neiDt_i(3)))^2)/121;
                    msc=msc+1;
                end
            end
            msrdTmp=msrdTmp(1:msc-1);
            
            
            if size(msrdTmp,1)>2
                msrdMn=mean(msrdTmp);
                msrdStd=std(msrdTmp);
                pdt_ind=pPc{ii}(jj,:);
    
                pdt_ind(1:3)=pdt_ind(1:3)-origin;
                pdt_ind(1:3)=pdt_ind(1:3)/len_scale;
                pdt_ind(2)=-pdt_ind(2);
    
                msrdDt(dsc,:)=[pdt_ind,...
                    1/2*(1+erf((1-msrdMn)/sqrt(2)/msrdStd))];
                dsc=dsc+1;
            end
            
        end
    end
    
    toc;
    
    msrdDt=msrdDt(1:dsc-1,:);
    
    dvBndDt=pdvb-origin;
    dvBndDt=dvBndDt/len_scale;
    dvBndDt(:,2)=-dvBndDt(:,2);
    
    edBndDt=pedb-origin;
    edBndDt=edBndDt/len_scale;
    edBndDt(:,2)=-edBndDt(:,2);
    
    pRfDt=pRf-origin;
    pRfDt=pRfDt/len_scale;
    pRfDt(:,2)=-pRfDt(:,2);
    
    
    odr=sprintf('%s%s/msrd_one_lensc_result/',idr,csNm);
    mkdir(odr);
    
    save(sprintf('%s%s_msrd_data.mat',odr,csNm),'msrdDt');
    
    save(sprintf('%s%s_msrd_dvb_data.mat',odr,csNm),'dvBndDt');
    save(sprintf('%s%s_msrd_edb_data.mat',odr,csNm),'edBndDt');
    save(sprintf('%s%s_msrd_ref_data.mat',odr,csNm),'pRfDt');

end