ssf time clipping

matlab/+SSF/drawKymosCondition.m

@@ -18,15 +18,22 @@
 % 
 % cbLabel = 'ERK SSF';
 % targetChannels = 1;
-function drawKymosCondition(flist,classList,labelNames,outfile,cbLabel,targetChannels,stride)
+function drawKymosCondition(flist,classList,labelNames,outfile,cbLabel,targetChannels,stride,framesPerHour)
 
+if ~exist('xticksLabels','var')
+    xticksLabels.xticks = [];
+    xticksLabels.xlabels = {};
+end
 xClasses = unique(classList);
 
 nKeep = 0;
 clipRange = [2.5,97.5];
+% clipRange = [0.5,99.5];
 clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep);
 
 xyim = []; % one per flist
+xticks = [];
+xtickLabels = {};
 for cc = 1:length(xClasses)
     idxCC = find(strcmp(classList,xClasses(cc)));
     flistCC = flist(idxCC);
@@ -74,29 +81,43 @@ for cc = 1:length(xClasses)
         plot([x,x],[y,y+szim(ff,1)],'-k')
         % don't draw to your right...
 %         plot([x+szim(ff,2),x+szim(ff,2)],[y,y+szim(ff,1)],'-k')
-        text(x,y,ccLabels{ff},'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.5],...
+        text(x,y,ccLabels{ff},'Interpreter','none','FontSize',6,'BackgroundColor',[0.9,0.5,0.7,0.5],...
             'VerticalAlignment','top');
+        
+        % label time
+        tHours = [24:24:szim(ff,2)-12];
+        for i = 1:length(tHours)
+            text(x+tHours(i),y+szim(ff,1),num2str(round(tHours(i)./framesPerHour)),'Interpreter','none','FontSize',6,...
+            'VerticalAlignment','bottom');
+        end
+        
+        
+
+
         x = x + szim(ff,2);
 %         date_movie = regexp(flistCC(ff).name,'(\d{4}-\d{2}-\d{2}).+_(\d{2})','tokens');
 %         text(x,y,[date_movie{1}{1} '_' date_movie{1}{2}],'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.5]);
 %         plot([x+62,x+62],[y+5,y+size(imq,1)-5],'-r')
 %         treatment = regexp(tblCC.Treatment{ff},'(.*)\w*(','tokens');
 %         text(x+30,y+20,[date_movie{1}{1} '_' date_movie{1}{2} ' :: ' treatment{1}{1}],'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.7]);
+        
     end
 
+
     % set background (index 0, colormap(1)) to white
     cm = parula();    
     cm(1,:) = [1,1,1];
     colormap(cm);
     axis on
-    set(gca,'xtick',[]);
+    set(gca,'xtick',xticksLabels.xticks);
+    set(gca,'XTickLabel',xticksLabels.xlabels);
     set(gca,'ytick',[]);
-    ylabel('xy \Downarrow x');
+    ylabel('xyz \Downarrow x');
     xlabel('time');
     fx = strfind(flist(ff).name,'.LEVER');
     fx = fx(1)-1;    
-    title(xClasses{cc},'Interpreter','none');
-    tick_labels = round(linspace(clipLimits(1),clipLimits(2),5),2);
+%     title(xClasses{cc},'Interpreter','none');
+    tick_labels = round(linspace(clipLimits(1),clipLimits(2),5),2,'significant');
     cb = colorbar('Ticks',[1,64,128,192,255],'TickLabels',tick_labels);
     cb.Label.String = cbLabel;
     if 1 == cc

matlab/+SSF/draw_ssf_ncd.m

@@ -2,23 +2,20 @@
 ROOT = '/g/leverjs/Olivier/Agne/march_2023_20x';
 flist = dir(fullfile(ROOT,'*.LEVER'));
 
-className = {};
-for ff=1:length(flist)
-    cn = regexp(flist(ff).name,'(\d+-\d+-\d+)_.*(\d\d)_ori*','tokens');
-    cn = [cn{1}{1} '_' cn{1}{2}];
-    className{ff} = classMap(cn);
-end
+classList = get_ssClasses(flist);
+[classUnique,ia,ic]= unique(classList);
+idxTrue = ic;
 % targetClass = {'AKT1_E17K','WT','PIK3CA_E545K'}
 % idx = find(cellfun(@(x) ~isempty(find(strcmp(targetClass,x))),className));
 % A=d(idx,idx);
 % className = className(idx);
 
 A = d;
-K=3;
-
+K=6;
+[idx,Y] = Cluster.SpectralCluster(A,K);
 
+[csf_mean, csf_std] = CSF.csf_spatial(Y,idxTrue)
 
-[idx,Y] = Cluster.SpectralCluster(A,K);
 % Y = tsne(Y,'Algorithm','exact','NumDimensions',3);
 clf;hold off;
 % if size(Y,2) == 3
@@ -26,21 +23,17 @@ clf;hold off;
 % else
 %     plot(Y(:,1),Y(:,2),'.')
 % end
-sym = {'ko','cv','m^','r*','gx','<k','b>','gs','bd'};
+sym = {'ko','cv','m^','r*','gx','gs','bd'};
 
 
-classList = unique(className);
-hx = zeros(size(classList));
+hx = zeros(size(classUnique));
 hold on
 
 for ff=1:size(Y,1)
     
-    idxSym = find(strcmp(classList,className{ff}));
-    if K==3
+    idxSym = find(strcmp(classUnique,classList{ff}));
+    
         h1 = plot3(Y(ff,1),Y(ff,2),Y(ff,3),sym{idxSym});
-    else
-        h1 = plot(Y(ff,1),Y(ff,2),sym{idxSym},'UserData',ff);
-    end
     if ~hx(idxSym)
         hx(idxSym) = h1;
     end
@@ -52,5 +45,8 @@ for ff=1:size(Y,1)
     %     end
 end
 set(gcf,'color','w')
-xlabel('u1');ylabel('u2'),zlabel('u3')
-legend(hx(find(hx)),classList(find(hx)),'interpreter','none');
\ No newline at end of file
+xlabel('k1');ylabel('k2'),zlabel('k3')
+legend(hx(find(hx)),classUnique(find(hx)),'interpreter','none');
+tstr = ['CSF^6_{Erk,Velocity} = ' mat2str([csf_mean, csf_std],2)]; 
+title(tstr)
+

matlab/+SSF/getClipLimits.m

 % map [-1,1] kymographs to color space [0,255]
 % 
-function clipLimits = getClipLimits(flist,targetChannelNumbers,clipRange,nKeep)
+function clipLimits = getClipLimits(flist,targetChannelNumbers,clipRange,nKeep,tClip)
 
 if ~exist('clipRange','var')
     clipRange= [2.5,97.5];
@@ -8,13 +8,16 @@ end
 if ~exist('nKeep','var')
     nKeep = 0;
 end
+if ~exist('tClip','var')
+    tClip = [];
+end
 p = ljsStartParallel(); 
 kymoPixels = {};
 parfor ff=1:length(flist)
     kp = {};
     for c = 1:length(targetChannelNumbers)
         strDB_ff = fullfile(flist(ff).folder,flist(ff).name);
-        im_ff = SSF.loadImage(strDB_ff,targetChannelNumbers(c));
+        im_ff = SSF.loadImage(strDB_ff,targetChannelNumbers(c), 0, tClip);
         if nKeep>0
             im_ff(im_ff<0)=0;
         end

matlab/+SSF/go_ssf_ncd.m

@@ -4,7 +4,7 @@
 % AKT_CHANNEL = 2;
 % H2B_CHANNEL = 3;
 % targetChannels = [ERK_CHANNEL,AKT_CHANNEL];
-function d=go_ssf_ncd(kROOT,targetChannels,clipRange, nKeep)
+function d=go_ssf_ncd(kROOT,targetChannels,clipRange, nKeep, tClip)
 
 tStart = tic();
 
@@ -15,7 +15,15 @@ end
 if ~exist('clipRange','var')
     clipRange = [2.5,97.5];
 end
+
+if ~exist('tClip','var')
+    tClip = [];
+end
+if isstring(kROOT) || ischar(kROOT)
     flist = dir(fullfile(kROOT,'*.LEVER'));
+else
+    flist = kROOT;
+end
 d = zeros(length(flist));
 kymoPixels = {};
  
@@ -24,7 +32,7 @@ p = ljsStartParallel();
 if iscell(clipRange)
     clipLimits = clipRange;
 else
-    clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep);
+    clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep,tClip);
 end
 %
 cmdList = NCD.dParallelCommandList(ones(length(flist)),p.NumWorkers);
@@ -44,7 +52,7 @@ parfor i=1:H
         end
         strDB_ff = fullfile(flist(ff).folder, flist(ff).name);
         if ~strcmp(strDB_ff,str_ff_prev)
-            im_ff = SSF.loadImage(strDB_ff, targetChannels, nKeep);
+            im_ff = SSF.loadImage(strDB_ff, targetChannels, nKeep, tClip);
             if iscell(clipLimits)
                 cx = clipLimits{ff};
             else 
@@ -59,7 +67,7 @@ parfor i=1:H
             dxx(i,j) = d1;
         else
             strDB_gg = fullfile(flist(gg).folder,flist(gg).name);
-            im_gg = SSF.loadImage(strDB_gg,targetChannels, nKeep);
+            im_gg = SSF.loadImage(strDB_gg,targetChannels, nKeep, tClip);
             if iscell(clipLimits)
                 cx = clipLimits{gg};
             else 
@@ -74,7 +82,7 @@ parfor i=1:H
             dxx(i,j) = min(d1,d2);
         end
         tElapsed = toc(t0);
-        fprintf('%s tElapsed = %0.2f\n',mat2str([ff,gg,dxx(i,j)]),tElapsed);
+%         fprintf('%s tElapsed = %0.2f\n',mat2str([ff,gg,dxx(i,j)]),tElapsed);
     end
 end
 % convert dxx to distance matrix d

matlab/+SSF/loadImage.m

 
-function im = loadImage(strDB, channelList, nKeep)
+function im = loadImage(strDB, channelList, nKeep, tClip)
 
 if ~exist('bDouble','var')
     bDouble = false;
@@ -9,6 +9,11 @@ time = 1; % kymo is single time point
 for c=1:length(channelList)
     im(:,:,:,c) = leversc.loadImage(strDB,time,channelList(c));  
 end
+% tClip
+if exist('tClip','var') && ~isempty(tClip)    
+    tClip(tClip>size(im,3)) = [];
+    im = im(:,:,tClip,:);
+end
 
 if ~exist('nKeep','var')
     return
@@ -20,3 +25,4 @@ if nKeep < 0
     im(im>0) = 0;
     im = abs(im);
 end
+

matlab/+Segment/frameSegment.m

@@ -143,7 +143,7 @@ for n=1:max(L(:))
     newCell.LoG = LoG.sampleCell(newCell,CONSTANTS,imLoG);
 
     % check ecc
-    v = eig(cov(newCell.pts));
+    v = eig(cov(double(newCell.pts)));
     ecc = max(v)/min(v);
     if ecc>50 
         continue