Merge branch 'log_seg' of git-bioimage.coe.drexel.edu:bioimage/leverjs into log_seg

leversc/src/MATLAB/@leversc/leversc.m

@@ -30,6 +30,7 @@ classdef leversc <handle
                 optional_figNumber=1;
             end
             [im,imD,figNumber]=parseArgs(optional_im,optional_imD,optional_figNumber);
+            im = imPad4(im);
             obj.figNumber=figNumber;
             if isfield(imD,'imageData')
                 CONSTANTS=imD;
@@ -122,3 +123,19 @@ end
 
 end
 
+function im = imPad4(im)
+
+pad = [0,0,0];
+for d = 1:3
+    if 0~=mod(size(im,d),4)
+        pad(d) = size(im,d) + 4-mod(size(im,d),4);
+    end
+end
+if any(pad)
+    sz = size(im);
+    sz = sz(1:3);
+    pad = max(pad,sz);
+    im(pad(1),pad(2),pad(3),:) = 0;
+end
+
+end
\ No newline at end of file

matlab/+Batch/batchProc.m

 
 function batchProc(ROOT,nproc,nthread)
 
-parpool(nproc)
+ljsStartParallel(nproc);
 
 flist = dir(fullfile(ROOT,'*.LEVER'));
 parfor ff=1:length(flist)

matlab/+Batch/processMovie.m

@@ -49,7 +49,8 @@ processInfo.tPatch=getElapsed(tx);
 
 tx=tic();
 % nMerged = Smooth.Merge.goMerge(conn,CONSTANTS,segParams,nProcessors);
-nMerged = Smooth.Merge.goMerge2(conn, CONSTANTS, segParams);
+% nMerged = Smooth.Merge.goMerge2(conn, CONSTANTS, segParams);
+nMerged = -1;
 processInfo.tMerge=getElapsed(tx);
 
 tx=tic();

matlab/+DenoiseCache/denoiseFolder.m

@@ -8,7 +8,7 @@ if 2==length(nProcessors)
     nProcessors=nProcessors(2);
 end
 
-nProcessors = min(nProcessors,16);
+nProcessors = min(nProcessors,8);
 
 flist=dir(fullfile(folderName,'*.h5'));
 if 0==length(flist)

matlab/+Helpers/cellEccentricity.m

+function ecc = cellEccentricity(cx)
+
+evals = eig(cov(cx.pts));
+ecc = min(evals)/max(evals);
\ No newline at end of file

matlab/+Lineage/eraseLineage.m

+function eraseLineage(conn)
+
+exec(conn,'delete from tblFamilies;delete from uiExtFamilies');
+strDB = conn.DataSource;
+Batch.batchTrack(strDB);
+Smooth.Patch.patch11(conn);
+Smooth.Patch.patch12(conn);
+Smooth.Patch.patch21(conn);
+Smooth.Patch.patchOcclusion(conn);
\ No newline at end of file

matlab/+MovieMaker/ssfMovie_erk.m

@@ -5,8 +5,6 @@ function ssfMovie1(strDB)
 bPreview = true
 
 
-
-
 bwx = bwLoG_pos;
 bwx(:,:,:,2) = bwLoG_neg;
 imd = MicroscopeData.MakeMetadataFromImage(bwx);

matlab/+Read/findImageFolder.m

@@ -42,7 +42,7 @@ if ~exist(ljs,'file')
 end
 [fid,message]=fopen(fullfile(leverFolder,'lever.json'));
 if (fid<=0)
-    fprintf(1,'ERROR -- getConstants.m -- image path invalid and no lever.json ROOT override (message=%s)\n',message);
+    fprintf(1,'ERROR %s -- getConstants.m -- image path invalid and no lever.json ROOT override (message=%s)\n',CONSTANTS.imageData.DatasetName,message);
     return;
 end
 

matlab/+Read/getCell.m

 function Cell=getCell(conn,cellID,CONSTANTS)
 
 Cell=[];
-
+qTable=fetch(conn,'pragma table_info(tblCells)');
+bLoG = any(cellfun(@length,(strfind(qTable.name,'LoG'))));
 try
-    cmd=['SELECT cellID,time,trackID,centroid,u16pixels,maxRadius,channel,area,surface,segCC '...
-        ' FROM tblCells WHERE cellID=' num2str(cellID)];
+    cmd=['SELECT cellID,time,trackID,centroid,u16pixels,maxRadius,channel,area,surface,segCC '];
+    if bLoG
+        cmd = [cmd ' ,LoG '];
+    end
+    cmd =[cmd ' FROM tblCells WHERE cellID=' num2str(cellID)];
     Q=ljsFetch(conn,cmd);
 catch
     try
         Cell=Read.getCell3D(conn,cellID);
     catch
         fprintf(2,'unable to read cell from %s\n',CONSTANTS.imageData.DatasetName);
-        return
     end
+    return
 end
 if isempty(Q)
     return
@@ -32,6 +36,9 @@ Cell.area=Q{8};
 Cell.surface=jsondecode(Q{9});
 Cell.segCC=Q{10};
 
+if bLoG
+    Cell.LoG = Q{11};
+end
 if exist('CONSTANTS','var')
     Cell = Read.setCellsIdxPts(Cell,CONSTANTS);
 end

matlab/+Read/getCell3D.m

@@ -2,7 +2,14 @@ function Cell=getCell3D(conn,cellID,CONSTANTS)
 
 Cell=[];
 
-cmd=['SELECT time,trackID,centroid,channel,area,maxRadius,verts,faces,edges,normals,u16pixels FROM tblCells WHERE cellID=' num2str(cellID)];
+qTable=fetch(conn,'pragma table_info(tblCells)');
+bLoG = any(cellfun(@length,(strfind(qTable.name,'LoG'))));
+
+cmd=['SELECT time,trackID,centroid,channel,area,maxRadius,verts,faces,edges,normals,u16pixels '];
+if bLoG
+    cmd = [cmd ',LoG '];
+end
+cmd = [cmd ' FROM tblCells WHERE cellID=' num2str(cellID)];
 Q=ljsFetch(conn,cmd);
 % cellid time trackID centroid channel area, maxRadius,
 % verts,edges,normals,faces,u16pixels

matlab/+SSF/catKymographs.m

+% pads time between images
+function i12 = catKymographs(i1,i2,dim)
+[i1,i2] = padIm(i1,i2);
+i12 = cat(dim,i1,i2);
+4;
+
+function [i1,i2] = padIm(i1,i2)
+
+if size(i1,3)>size(i2,3)
+    i2(:,:,end:size(i1,3),:)=0;
+else
+    i1(:,:,end:size(i2,3),:)=0;
+end
+
+if size(i1,2)>size(i2,2)
+    i2(:,end:size(i1,2),:,:)=0;
+else
+    i1(:,end:size(i2,2),:,:)=0;
+end
+
+if size(i1,1)>size(i2,1)
+    i2(end:size(i1,1),:,:,:)=0;
+else
+    i1(end:size(i2,1),:,:,:)=0;
+end
+

matlab/+SSF/cropKymo.m

+function im_crop = cropKymo(imSrc,idxPixels)
+
+% first copy over target pixels
+im = 0 * imSrc;
+im(idxPixels) = imSrc(idxPixels);
+
+[r,c,z]=ind2sub(size(im),idxPixels);
+
+r = cropDim(r,size(im,1));
+c = cropDim(c,size(im,2));
+z = cropDim(z,size(im,3));
+
+im_crop = im(min(r):max(r),min(c):max(c),min(z):max(z),:);
+
+
+function y = cropDim(x,sz)
+
+y = [min(x)-1:max(x)+1];
+y = min(y,sz);
+y = max(y,1);
+
+

matlab/+SSF/drawMitoticSSFtracks.m

+% erk1D = Y(:,1);
+% [~,idxSort] = sort(erk1D);
+% qMitoSorted = qMitotic(idxSort,:);
+% 
+function imx = drawMitoticSSFtracks(qMitotic,erkClipLimits,bDraw)
+
+if ~exist('bDraw','var')
+    bDraw = false;
+end
+
+for i = 1:height(qMitotic)
+    if ~isa(qMitotic.erk{i},'uint8')
+        qMitotic.erk{i} = SSF.quantize8(qMitotic.erk{i},erkClipLimits);
+    end
+end
+
+tmax = max(cellfun(@length,qMitotic.erk));
+imx = zeros(2*height(qMitotic),tmax);
+for i = 1:height(qMitotic)
+    erk_i = qMitotic.erk{i};
+    imx(2*i,end-length(erk_i)+1:end) = erk_i;
+end
+
+if ~bDraw
+    return
+end
+
+
+cmPos = colormap('parula');
+cmNeg = colormap('cool');
+cm = [cmNeg(1:2:end,:);cmPos(1:2:end,:)];
+cm(1,:) = [1,1,1];
+
+
+figure;imagesc(imx);colormap(cm);
+c = colorbar();
+clim = [min(c.Ticks),max(c.Ticks)];
+cticks = [clim(1), mean([clim(1),mean(clim)]), mean(clim), mean([clim(2),mean(clim)]),clim(2)];
+c.Ticks = cticks;
+c.TickLabels = {'BG','-0.5','0','0.5','1'}
+c.Label.String = 'Erk SSF';
+
+% close all
+
+% cm = colormap('parula');
+% cm(1:13,:)=repmat([1,1,1],13,1);
+% % c = colorbar('ticks',round(linspace(double(erkLimits(1)),double(erkLimits(2)),10)));
+% c = colorbar
+% c.Label.String = 'Erk SSF (8 bit, clipped to 95%)'
+% 
+% xx = round(linspace(50,size(imx,2),5));
+% set(gca,'XTickLabel',arrayfun(@(y) num2str((y-308),xx,'UniformOutput',false))
+% xlabel('frames before mitosis');
+% ylabel('mitotic parent track sorted by NCD principal component')
+% yy = yticks();
+% set(gca,'YTickLabel',arrayfun(@(y) num2str(y/2),yticks()))
+% 4;

matlab/+SSF/drawSSFclusters.m

+% erk1D = Y(:,1);
+% [~,idxSort] = sort(erk1D);
+% qMitoSorted = qMitotic(idxSort,:);
+% 
+% function imx = drawSSFtracks(qMitotic,idx,bDraw)
+
+% if ~exist('bDraw','var')
+%     bDraw = false;
+% end
+tmax = max(cellfun(@length,qMitotic.erk));
+imx = [];
+imxH = 1;
+clusterPts = [];
+for k = 1:max(idx)
+    idxk = find(idx==k);
+    qk = qMitotic(idxk,:);
+    qk.yk = Y(idxk,1);
+    qk = sortrows(qk,'yk');
+    for i = 1:height(qk)
+        erk_i = qk.erk{i};
+        imx(imxH,tmax-length(erk_i)+1:tmax) = erk_i;
+        imxH = imxH+2;
+    end
+    imxH = imxH + 2;
+    clusterPts=[clusterPts,imxH];
+    imxH = imxH + 4;
+
+end
+
+
+
+% if ~bDraw
+%     return
+% end
+
+
+close all
+
+cm = colormap('parula');
+cm(1,:) = [1,1,1];
+
+figure;imagesc(imx);colormap(cm);hold on
+c = colorbar('ticks',[2,64,128,192,255], 'ticklabels',{'-1.0','-0.5','0','0.5','1'});
+c.Label.String = 'Erk SSF';
+for p = 1:length(clusterPts)
+    plot([1,size(imx,2)],[clusterPts(p),clusterPts(p)],'-..k')
+end
+
+% 
+
+% cm = colormap('parula');
+% cm(1:13,:)=repmat([1,1,1],13,1);
+% % c = colorbar('ticks',round(linspace(double(erkLimits(1)),double(erkLimits(2)),10)));
+% c = colorbar
+% c.Label.String = 'Erk SSF (8 bit, clipped to 95%)'
+% 
+xx = round(linspace(50,size(imx,2),5));
+set(gca,'XTick',xx);
+set(gca,'XTickLabel',arrayfun(@(y) num2str(y-308),xx,'UniformOutput',false))
+xlabel('frames before mitosis');
+ylabel('mitotic parent track sorted by cluster and NCD principal component')
+yy = yticks();
+set(gca,'YTickLabel',arrayfun(@(y) num2str(y/2),yticks(),'uniformOutput',false))
+% h=gca; h.YAxis.TickLength = [0 0];
+yticks([])
+% 4;
+set(gcf,'color','w')
+% exportgraphics(gcf,'erk ssf heat map tracks NCD sorted cluster FLIF1d.pdf','Resolution',600,'contenttype','image')

matlab/+SSF/drawSSFtracks.m

+% strDB = '/g/leverjs/Olivier/Agne/march_2023_20x/2022-04-05_steady_state_48h_deprived_13_ori.LEVER';
+
+function drawSSFtracks(ROOT,cache, kymo_channel, channel_description,fnames)
+
+[~,dsName,ext] = fileparts(ROOT);
+if ~exist('fnames','var')
+    fnames = [];
+end
+
+outfile = [dsName ' ' channel_description ' SSF projection_y_t.pdf'];
+if isempty(ext)
+    flist = dir(fullfile(ROOT,'*.LEVER'));
+else
+    flist = dir(ROOT);
+end
+f=figure(1);
+f.WindowState = 'maximized';
+pause(1)
+
+kymoPixels = {};
+for ff=1:length(flist)
+    strDB = fullfile(flist(ff).folder,flist(ff).name);
+    [ROOT,lf] = fileparts(strDB);
+    strKymoDB = fullfile(ROOT,cache,[lf '.LEVER_ssf_cache.LEVER']);
+    if ~exist(strKymoDB,'file')
+        continue
+    end
+    imKymo = SSF.loadImage(strKymoDB,kymo_channel);
+    kymoPixels{ff} = imKymo(find(imKymo));
+end 
+kymoPixels = vertcat(kymoPixels{:});
+erkClipLimits = [prctile(kymoPixels,2.5),prctile(kymoPixels,97.5)];
+
+for ff = 1:length(flist)
+    strDB = fullfile(flist(ff).folder,flist(ff).name);
+    [ROOT,lf] = fileparts(strDB);
+    strKymoDB = fullfile(ROOT,cache,[lf '.LEVER_ssf_cache.LEVER']);
+    if ~exist(strKymoDB,'file')
+        continue
+    end
+    imKymo = SSF.loadImage(strKymoDB,kymo_channel);
+    % 
+    imp = squeeze(max(imKymo,[],2));
+    imn = squeeze(min(imKymo,[],2));
+    imn = abs(imn);
+    idx = find(imn>imp);
+    imp(idx) = -1 .* imn(idx);
+%     imp(imp<-0.5) = 0;
+    impq = SSF.quantize8(imp,erkClipLimits);
+%     impq(imp==0 & imn==0) = 0;
+    
+%     % keep just [0,1] SSF
+%     imp = squeeze(max(imKymo,[],2));
+%     erkClipLimits = [0,prctile(kymoPixels,99.9)];
+%     impq = SSF.quantize8(imp,erkClipLimits);
+%     impq(imp==0) = 0; % set the background
+%     
+    
+    clf;imagesc(impq)
+    cm = colormap('parula');
+    cm(1,:) = [1,1,1]; % make the background white!
+    colormap(cm)
+    c = colorbar();
+    clim = [min(impq(:)),max(impq(:))];
+    cticks = [clim(1)+1, mean([clim(1),mean(clim)]), mean(clim), mean([clim(2),mean(clim)]),clim(2)];
+    c.Ticks = cticks;
+    % c.TickLabels = {'-1','-0.5','0','0.5','1'}
+
+    ticksSSF = linspace(erkClipLimits(1),erkClipLimits(2),5);
+    ticksSSF = arrayfun(@(x) num2str(x,2),ticksSSF,'UniformOutput',false);
+    c.TickLabels = ticksSSF;
+    c.Label.String = [channel_description ' SSF (dataset reference)'];
+    if ~isempty(fnames)
+        lf = fnames{ff};
+    end
+    hold on
+    yl = ylim();
+    plot([62,62],[yl(1),yl(2)],'--r','linewidth',2)
+    title([lf ' ' channel_description ' SSF projection to (Y/2,time)'],'interpreter','none')
+    xlabel('time (frames)')
+    ylabel('Y (max intensity projection along X)')
+    bAppend = ff>1;
+    exportgraphics(gcf,outfile,'Resolution',600,'contenttype','image','append',bAppend);
+    
+end
\ No newline at end of file

matlab/+SSF/draw_ssf_ncd.m

 % load('agne_72_8bpp.mat')
 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
+% 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
 % targetClass = {'AKT1_E17K','WT','PIK3CA_E545K'}
 % idx = find(cellfun(@(x) ~isempty(find(strcmp(targetClass,x))),className));
 % A=d(idx,idx);
@@ -28,7 +28,7 @@ clf;hold off;
 sym = {'ko','cv','m^','r*','gx','<k','b>','gs','bd'};
 
 
-classList = unique(classMap.values);
+classList = unique(className);
 hx = zeros(size(classList));
 hold on
 

matlab/+SSF/draw_ssf_tracks.m

-close all
-classList = unique(classMap.values);
-for ff = 1:length(flist)
-    strDB = fullfile(flist(ff).folder,flist(ff).name);
-    [conn,CONSTANTS,segParams]=openDB(strDB);
-
-    cmd = 'select trackID from tblCells inner join tblFamilies on cellID = cellID_parent';
-    q = fetch(conn,cmd);
-    close(conn);
-
-    erk_tracks = movie_tracks{ff};
-    class = classMap(flist(ff).name(1:13));
-    id = find(cellfun(@length,strfind(classList,class)));
-    figure(id) ; hold on
-
-    tids = erk_tracks.keys;
-    tids = [tids{:}];
-    tids = intersect(tids,q.trackID);
-    for i = 1:length(tids)        
-        ssf = erk_tracks(tids(i));
-        if size(ssf,2) <= 10 || all(0==ssf(1,:))
-            continue
-        end
-        plot(5*ssf(2,:),ssf(1,:),'-');
-    end
-end
-
-for i = 1:length(classList)
-    figure(i)
-    set(gcf,'color','w')
-    title(classList{i},'Interpreter','none')
-    xlabel('time (minutes)')
-    ylabel('SSF')
-    fname = ['erk_tracks_' classList{i} '.tif'];
-    print(fname,'-dtiffn');
-end
-
-cx={}
-for ff = 1:length(flist)
-        cx{ff} = classMap(flist(ff).name(1:13))
-end
\ No newline at end of file

matlab/+SSF/gen_SSF_kymo.m

 
+% BE CAREFUL HERE! BE DRAGONS! GET THESE target names & numbers RIGHT OR PERISH
+% targetChannelNames = {'Erk SSF'; 'Akt SSF' ; 'DHB SSF'; 'H2B SSF'; 'velocity SSF'}; % channel names in the kymo
+% targetChannelNames = {'Erk SSF'; 'Oct4 SSF' ;  'H2B SSF'; 'velocity SSF'}; % channel names in the kymo
+% targetChannelNumbers = [2, 3, -1]; % maps image channel names to kymo channel names
+function gen_SSF_kymo(strDB,targetChannelNames,targetChannelNumbers,outfolder,bDownscale, startFrame)
 
-function gen_SSF_kymo(ROOT)
-targetChannelNames = {'Erk SSF'; 'Akt SSF' ; 'DHB SSF'; 'H2B SSF'; 'velocity SSF'}; % channel names in the kymo
-targetChannelNumbers = [-2, -3, 1, -4]; % maps image channel names to kymo channel names
+if ~exist('bDownscale','var')
+    bDownscale = false;
+end
+if ~exist('startFrame','var')
+    startFrame = 1;
+end
 
-cacheFolder = fullfile(ROOT,'kymoV2');
-if ~exist(cacheFolder,'dir')
-    mkdir(cacheFolder);
+[fx,nx,ext]=fileparts(strDB);
+if strcmp(ext,'.LEVER')
+    flist = dir(strDB);
+    ROOT = fx;
+else
+    ROOT = strDB;
+    flist = dir(fullfile(strDB,'*.LEVER'));
 end
-flist = dir(fullfile(ROOT,'*.LEVER'));
 L = length(flist);
 
+
+if ~exist(outfolder,'dir')
+    mkdir(outfolder);
+end
+
 imd={};
+p = ljsStartParallel(16);
+
 parfor ff=1:L
+% for ff=1:L    
     str_ff = fullfile(flist(ff).folder,flist(ff).name);     
-    if exist(fullfile(cacheFolder,[flist(ff).name '_ssf_cache.json']),'file')
+    if exist(fullfile(outfolder,[flist(ff).name '_ssf_cache.json']),'file')
         fprintf(1,'genKymo skipping %s kymo exists\n',str_ff);
         continue
     end
     t0 = tic();
-    imVolume = SSF.ssf_volume(str_ff,targetChannelNumbers);
+    imVolume = SSF.ssf_volume(str_ff,targetChannelNumbers, bDownscale);
     if isempty(imVolume)
         continue
     end
-    imVolume(:,:,:,5) = SSF.velocity_kymo(str_ff);
-    imVolume = imVolume(:,:,31:end,:); % trim off first 30 frames
+%     imVolume(:,:,:,length(targetChannelNumbers)+1) = SSF.velocity_kymo(str_ff);
+    imVolume = imVolume(:,:,startFrame:end,:); % trim off first 30 frames
     imd{ff}=MicroscopeData.MakeMetadataFromImage(imVolume);
     imd{ff}.DatasetName=[flist(ff).name '_ssf_cache'];
-    imd{ff}.imageDir = cacheFolder;
+    imd{ff}.imageDir = outfolder;
     imd{ff}.ChannelNames = targetChannelNames;
-    MicroscopeData.WriterH5(imVolume,'imageData',imd{ff},'path',cacheFolder);    
+    MicroscopeData.WriterH5(imVolume,'imageData',imd{ff},'path',outfolder);  
+
     tElapsed = toc(t0);
     fprintf(1,'ssf_volume %d tElapsed=%0.2f\n',ff,tElapsed);
 end
-Import.leverImport('',cacheFolder);
+Import.leverImport('',outfolder);
\ No newline at end of file

matlab/+SSF/goErkCN.m

+% strDB = '/f/leverjs/Olivier/Yannick/20210924/20210924_pos3.LEVER';
+% tblCN = goErkCN(strDB,174279)
+% tblCN = goErkCN(strDB,224589)
+function tblCN = goErkCN(strDB,targetTrack)
+
+[conn,CONSTANTS,segParams]=openDB(strDB);
+[ROOT,lf] = fileparts(strDB);
+strKymo = fullfile(ROOT,'kymoV',[lf '.LEVER_ssf_cache.LEVER']);
+imKymo = leversc.loadImage(strKymo,1,1);
+
+cmd = ['select cellID,trackID,centroid,time,LoG from tblCells where trackID = '...
+    num2str(targetTrack)];
+q=fetch(conn,cmd);
+
+erkSSF = [];
+for i = 1:height(q)
+    centroid = jsondecode(q.centroid{i});
+    centroid = [round(centroid(1:2))', q.time(i)];
+    sample_value = SSF.ssf_sample(imKymo,centroid);
+    erkSSF(i) = sample_value;
+end
+
+featureParams.channel = 2;
+featureParams.name='CytoToNuclearRatio';
+featureParams.color="#00ffff";
+tblCN = table();
+for i = 1:height(q)
+    t = q.time(i);
+    cellList = Read.getCellsTime(conn,t,segParams.channels(1),CONSTANTS);
+    nt = SSF.cytoToNuclearRatio(cellList,q.cellID(i),t,CONSTANTS,featureParams,segParams);
+    nt.erkSSF = erkSSF(i);
+    tblCN = [tblCN;nt]
+end
+ee = imKymo(find(imKymo));
+erkClipLimits = [prctile(ee,2.5), prctile(ee,97.5)];
+
+erk = tblCN.erkSSF;
+erk = SSF.quantize8(erk,erkClipLimits);
+erk = 2 * (double(erk) / 255) -1;
+figure;plot(erk)
+hold on
+plot(tblCN.cnRatio)
\ No newline at end of file

matlab/+SSF/goTrackSSFkymo.m

+
+ROOT = '/f/leverjs/Olivier/Yannick/20210924/';
+flist = dir(fullfile(ROOT,'*.LEVER'));
+
+cxKymo = {};
+qMitotic = {};
+
+for ff = 1:length(flist)
+    strDB = fullfile(flist(ff).folder,flist(ff).name);
+    [ROOT,lf] = fileparts(strDB);
+    strKymoDB = fullfile(ROOT,'kymoV',[lf '.LEVER_ssf_cache.LEVER']);
+    cxKymo{ff} = SSF.loadImage(strKymoDB,1);
+    [conn,CONSTANTS,segParams]=openDB(strDB);    
+    imErk = cxKymo{ff}(:,:,:,1);
+    qm = fetch(conn,'select * from (select cellID,trackID,count(time) as nt from tblCells where trackID in (select trackID from tblCells inner join tblFamilies on cellID_parent=cellID) group by trackID) where nt>23');
+    for i = 1:length(qm.trackID)
+        qcx = fetch(conn,['select centroid,time from tblCells where trackID = ' num2str(qm.trackID(i))]);
+        centroids = cellfun(@jsondecode,qcx.centroid,'UniformOutput',false);
+        centroids = reshape(vertcat(centroids{:}),3,[])';
+        centroids = round(centroids);  
+        qm.centroids{i} = [centroids(:,1:2),qcx.time];
+        qm.idxCentroids{i} = sub2ind(size(cxKymo{ff}),centroids(:,2),centroids(:,1),qcx.time);        
+        qm.movieID(i) = ff;
+        qm.erk{i} = SSF.ssf_sample_idx(imErk,qm.idxCentroids{i});
+        qm.isEmpty(i) = all(0==qm.erk{i});       
+    end
+    qm(qm.isEmpty,:)=[];
+    qMitotic = [qMitotic;qm];
+end
+% qMitotic.erk is the signal. quantize it to 8 bit now
+erkCombined =[qMitotic.erk{:}]' ;
+erkClipLimits = [prctile(erkCombined,2.5), prctile(erkCombined,97.5)];
+for c = 1:length(cxKymo)
+    cxKymo{c} = SSF.quantize8(cxKymo{c},erkClipLimits);
+end
+% % a bit slow, but steady. find distancefrom each mitotis
+for i = 1:height(qMitotic)
+
+    strDB = fullfile(flist(qMitotic.movieID(i)).folder,flist(qMitotic.movieID(i)).name);
+    [conn,CONSTANTS,segParams]=openDB(strDB);
+    q1 = fetch(conn,['select centroid,time from tblCells where cellID = ' num2str(qMitotic.cellID(i))]);
+    colony = Read.getCellsTime(conn,q1.time,4,CONSTANTS);
+    for c = 1:length(colony)
+        bw = zeros(CONSTANTS.imageData.Dimensions(2),CONSTANTS.imageData.Dimensions(1));
+        bw(colony(c).idxPts) = 1;
+        centroid1 = round(jsondecode(q1.centroid{1}));
+        if ~bw(centroid1(2),centroid1(1))
+            continue
+        end
+        bwd = bwdist(~bw);
+        bwd = bwd./max(bwd(:));
+        qMitotic.dNormalized(i) = bwd(centroid1(2),centroid1(1));
+        qMitotic.time(i) = q1.time;
+    end
+end
+
+
+
+for i = 1:height(qMitotic)    
+    erk_i = medfilt1(qMitotic.erk{i}, 5);
+    erk_i = SSF.quantize8(erk_i, erkClipLimits);
+    qMitotic.erk{i} = erk_i;
+end
+
+d = [];
+L = height(qMitotic);
+parfor i = 1:L
+%     im_i = SSF.cropKymo(cxKymo{qMitotic.movieID(i)},qMitotic.idxCentroids{i});
+    for j = 1:L
+%         im_j = SSF.cropKymo(cxKymo{qMitotic.movieID(j)},qMitotic.idxCentroids{j});
+        erk_j = medfilt1(qMitotic.erk{j},5)
+        im_j = SSF.quantize8(erk_j,erkClipLimits);
+        d(i,j) = SSF.ncd_ssf_volume(im_i,im_j);
+    end
+end
+
+% for i = 1:length(d)
+%     for j = 1 : length(d)
+%         d(i,j)=max(d(i,j),d(j,i));
+%     end
+%     d(i,i) = max(d(i,i),0);
+% end
+
+close all
+A = d;
+K=3;
+[idx,Y] = Cluster.SpectralCluster(A,K);
+% Y(:,3) = [erkMitotic.d];
+figure;hold on
+hred = [];
+hgreen = [];
+for i = 1 : height(qMitotic)    
+    if qMitotic.dNormalized(i)<.2929
+        hred = [hred,plot3(Y(i,1),Y(i,2),Y(i,3),'.r')];        
+    else
+        hgreen = [hgreen,plot3(Y(i,1),Y(i,2),Y(i,3),'.g')];        
+    end
+end
+xlabel('u1')
+ylabel('u2')
+zlabel('u3')
+legend([hred(1),hgreen(1)],'d_{normalized} < 0.29','d_{normalized} \geq .2929')
+
+set(gcf,'color','w')
+figure;hold on
+for i = 1 : height(qMitotic)
+    if qMitotic.dNormalized(i)<.2929
+        plot(Y(i,1),Y(i,2),'.r');
+    else
+        plot(Y(i,1),Y(i,2),'.g');
+    end
+end
+xlabel('u1')
+ylabel('u2')
+legend({'d_{normalized}>.2929','d_{normalized}<.2929'})
+set(gcf,'color','w')
+% exportgraphics(gcf,'erkMitoticParentTracksEmbedding_medfilt.pdf','Resolution',300,'contenttype','image')
\ No newline at end of file