Merge branch 'master' of https://git-bioimage.coe.drexel.edu/bioimage/leverjs

matlab/+Batch/batchFeatures.m

+function batchFeatures(conn,CONSTANTS,nProcessors)
+
+algorithms=Read.getAlgorithms(conn,'feature');
+if isempty(algorithms)
+    return
+end
+segParams=Read.getSegmentationParams(conn);
+CONSTANTS=CellFeatures.resetFeatures(conn,CONSTANTS,algorithms);
+
+if ~exist('nProcessors','var') || isempty(nProcessors)
+    nProcessors=feature('numcores');
+end    
+p=gcp('nocreate');
+if isempty(p) || p.NumWorkers~=nProcessors
+    delete(p);
+    p=[];
+end
+if isempty(p)
+    p=parpool(nProcessors);
+end
+
+strDB=conn.DataSource;    
+for iAlgorithm=1:length(algorithms)
+    tNext=1;
+    for idp=1:ceil(CONSTANTS.imageData.NumberOfFrames/p.NumWorkers)
+        featureTable=Composite();
+        spmd
+            tTarget=tNext+labindex-1; % labIndex is 1 based
+            if tTarget<=CONSTANTS.imageData.NumberOfFrames
+                cx=openDB(strDB);
+                featureTable=algorithms(iAlgorithm).function(cx,tTarget,CONSTANTS,algorithms(iAlgorithm).params,segParams);
+                close(cx);
+            else
+                featureTable=[];
+            end
+        end
+        for i=1:length(featureTable)
+            if ~isempty(featureTable{i})
+                CellFeatures.writeFeatureTable(conn,featureTable{i});
+            end
+        end
+        tNext=tNext+p.NumWorkers;
+    end
+
+end
+
+4;
\ No newline at end of file

matlab/+Batch/processMovie.m

@@ -22,5 +22,8 @@ tx=tic();
 Smooth.Classify.goClassifyMitosis(conn,CONSTANTS);
 Smooth.extFamily(conn);
 tProcess.tClassify=toc(tx);
+tx=tic();
+Batch.batchFeatures(conn,CONSTANTS);
+tProcess.tFeatures=toc();
 
 close(conn)

matlab/+CellFeatures/cytoToNuclearRatio.m

-function cytoToNuclearRatio(conn,CONSTANTS,args,t)
+function tblCellFeatures=cytoToNuclearRatio(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
 
 params.channel=1;
 params.name='CytoToNuclearRatio';
 params.color="#00ffff";
 % check for param override
-params=Segment.getParams(params,args);
-
-% segmentation params -- for min radius
-algs=Read.getAlgorithms(conn,'segment');
-args=algs.params;
-segParams=Segment.getDefaultSegParams();
-segParams=Segment.getParams(segParams,args);
-
-rgRadius=Ensemble.getEnsembleRadius(segParams.minimumRadius_um);
-resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
-min_radius_pixels = rgRadius ./ resolution_um(1);
+featureParams=Segment.getParams(params,featureParams);
 
-% todo -- handle anistropy in 3d!
-
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
-if isempty(cellIDs)
-    return
+cellFeatures=[];
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+    return;
 end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
+
 im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
-    params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+    featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
 if all(isnan(im(:)))
     return;
 end
 is3d=size(im,3)>1;
 labelCells=0*im;
 dInterior=labelCells;
-for i=1:length(cellIDs)
-    if is3d
-        c=Read.getCell3D(conn,cellIDs(i));
-        idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1),c.pts(:,3));
-    else
-        c=Read.getCell(conn,cellIDs(i));
-        idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
-    end
+for i=1:length(cellList)
+    idx=cellList(i).idxPts;
     bwInterior=logical(0*im);
     bwInterior(idx)=1;
+    % we need bwdist for just this cell
     bwd=bwdist(~bwInterior);
+    % isolated bwdist put in combo image, indexed by labelCells
     dInterior=dInterior+bwd;
-    labelCells(idx)=c.cellID;
+    labelCells(idx)=cellList(i).cellID;
 end
 
 % find region surrounding segmentation
 [dExterior,idxd]=bwdist(labelCells);
+% 
+for i=1:length(cellList)
     
-for i=1:length(rgRadius)
-    li=labelCells;
-    lo=labelCells(idxd);
-    radius=rgRadius(i);
-    lo(dExterior>radius)=0;
-    % should this be 0, or min_radius_pixels? i.e. buffer between?
-    lo(dInterior>=radius)=0;
-    li(dInterior<radius)=0;
-    
-    labelOuter(:,:,:,i)=lo;
-    labelInner(:,:,:,i)=li;
-end
-
-for i=1:length(cellIDs)
-    if is3d
-        c=Read.getCell3D(conn,cellIDs(i));
-    else
-        c=Read.getCell(conn,cellIDs(i));
-    end
-    if isfield(c,'segCC')
-    iRadius=c.segCC-floor(c.segCC);
-        if 0==iRadius
-            iRadius=1;
-        else
-            prec=ceil(log10(length(rgRadius)+1));% # of digits in the radius tag
-            iRadius=iRadius*10^prec;
-            iRadius=int32(iRadius);
-        end
-    else
-        iRadius=1;
-    end
-    
-    idxInner=find(labelInner(:,:,:,iRadius)==c.cellID);
-    idxOuter=find(labelOuter(:,:,:,iRadius)==c.cellID);
+    % interior pixels > 2 pixels from cell surface
+    idxInner = find( labelCells==cellList(i).cellID & dInterior>=2);
+    % pixels <=2 pixels from cell surface 
+    idxInteriorBoundary=find( labelCells==cellList(i).cellID & dInterior<2);    
+    %  exterior pixels nearest to this cell 
+    idxExterior=idxd(find( labelCells(idxd)==cellList(i).cellID & dExterior<=3 & dExterior>0));
+    idxOuter=[idxInteriorBoundary;idxExterior];
     ratio=getErkRatio(im,idxOuter,idxInner);
-    if isnan(ratio)
-        cmd=['delete from tblCellFeatures where cellID=' num2str(cellIDs(i))];
-        exec(conn,cmd);
-    else
-        ratio=round(ratio,4);
-        cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
-            num2str(cellIDs(i)) ')'];
-        exec(conn,cmd);
-        cmd=['update tblCellFeatures set ' params.name '=' num2str(ratio) ...
-            ' where cellID=' num2str(cellIDs(i))];    
-        exec(conn,cmd);
-    end
-    
+    cellFeatures=[cellFeatures;cellList(i).cellID,ratio]; 
 end
-
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;
 
 function ratio=getErkRatio(im,idxOuter,idxInner)
 ratio=NaN;

matlab/+CellFeatures/netIntensity.m

-function netIntensity(conn,CONSTANTS,args,t)
+
+function tblCellFeatures=netIntensity(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
+cellFeatures=[];
 
 params.channel=1;
 params.name='NetIntensity';
 params.color="#ff0000";
 % check for param override
-params=Segment.getParams(params,args);
+featureParams=Segment.getParams(params,featureParams);
 
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+    return;
+end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
 
 im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
-    params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+    featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
 if all(im(:)==0)
     return
 end
 bias=median(im(:));
-for i=1:length(cellIDs)
-    c=Read.getCell(conn,cellIDs(i));
-    idxSegmentation=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
-    
-    net=sum(im(idxSegmentation)-bias);
-    cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
-        num2str(cellIDs(i)) ')'];
-    exec(conn,cmd);
-    cmd=['update tblCellFeatures set ' params.name '=' num2str(net) ...
-        ' where cellID=' num2str(cellIDs(i))];    
-    exec(conn,cmd);
+for i=1:length(cellList)    
+    net=sum(im(cellList(i).idxPts)-bias);
+    cellFeatures=[cellFeatures;cellList(i).cellID,net];
 end
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;

matlab/+CellFeatures/netNuclearEnvelope.m

-function netNuclearEnvelope(conn,CONSTANTS,args,t)
+
+function tblCellFeatures=netNuclearEnvelope(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
 
 params.channel=1;
 params.name='netNuclearEnvelope';
 params.color="#ff0000";
 % check for param override
-params=Segment.getParams(params,args);
-
-% segmentation params -- for min radius
-algs=Read.getAlgorithms(conn,'segment');
-args=algs.params;
-segParams=Segment.getDefaultSegParams();
-segParams=Segment.getParams(segParams,args);
-
-rgRadius=Ensemble.getEnsembleRadius(segParams.minimumRadius_um);
-resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
-min_radius_pixels = rgRadius ./ resolution_um(1);
-
-% todo -- handle anistropy in 3d!
+featureParams=Segment.getParams(params,featureParams);
 
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+    return;
+end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
 
 im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
-    params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+    featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
 if any(isnan(im))
     return;
 end
-is3d=size(im,3)>1;
 labelCells=0*im;
 dInterior=labelCells;
+% get separate distance transform for each cell
 for i=1:length(cellIDs)    
-    if is3d
-        c=Read.getCell3D(conn,cellIDs(i));
-        idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1),c.pts(:,3));
-    else
-        c=Read.getCell(conn,cellIDs(i));
-        idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
-    end
     bwInterior=logical(0*im);
-    bwInterior(idx)=1;
+    bwInterior(cellList(i).idxPts)=1;
     bwd=bwdist(~bwInterior);
     dInterior=dInterior+bwd;
-    labelCells(idx)=c.cellID;
+    labelCells(idx)=cellList(i).cellID;
 end
 
 % find region surrounding segmentation
 [dExterior,idxd]=bwdist(labelCells);
 
-for i=1:length(rgRadius)
-    li=labelCells;
-    lo=labelCells(idxd);
-    radius=rgRadius(i);
-    lo(dExterior>radius)=0;
-    % should this be 0, or min_radius_pixels? i.e. buffer between?
-    lo(dInterior>=radius)=0;
-    li(dInterior>=radius)=0;
-    
-    labelOuter(:,:,:,i)=lo;
-    labelInner(:,:,:,i)=li;
-end
 bias=median(im(:));
 for i=1:length(cellIDs)
-    if is3d
-        c=Read.getCell3D(conn,cellIDs(i));
-    else
-        c=Read.getCell(conn,cellIDs(i));
-    end
-    if isfield(c,'segCC')
-    iRadius=c.segCC-floor(c.segCC);
-        if 0==iRadius
-            iRadius=1;
-        else
-            prec=ceil(log10(length(rgRadius)+1));% # of digits in the radius tag
-            iRadius=iRadius*10^prec;
-            iRadius=int32(iRadius);
-        end
-    else
-        iRadius=1;
-    end
 
-    idxInner=find(labelInner(:,:,:,iRadius)==c.cellID);
-    idxOuter=find(labelOuter(:,:,:,iRadius)==c.cellID);
+    % read signal within 2 pixels on either side of the boundary
+    idxInner = find( labelCells==cellList(i).cellID & dInterior<=2);
+    idxOuter = idxd(find( labelCells(idxd)==cellList(i).cellID & dExterior<=2 & dExterior>0));   
     net=mean( [im(idxInner)-bias;im(idxOuter)-bias]);    
-    
     net=round(net,4);
-    cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
-        num2str(cellIDs(i)) ')'];
-    exec(conn,cmd);
-    cmd=['update tblCellFeatures set ' params.name '=' num2str(net) ...
-        ' where cellID=' num2str(cellIDs(i))];
-    exec(conn,cmd);    
-end
-
-
-function ratio=getErkRatio(im,idxOuter,idxInner)
-ratio=NaN;
-imOuter=im(idxOuter);
-imOuter(0==imOuter)=[];
-if length(imOuter)<5
-    return
-end
-imInner=im(idxInner);
-imInner(0==imInner)=[];
-if length(imInner)/length(idxInner)<0.1
-    return
+    cellFeatures=[cellFeatures;cellList(i).cellID,net];
 end
-ratio=median(imOuter)/median(imInner);
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;

matlab/+CellFeatures/resegFeatures.m

@@ -24,6 +24,7 @@ else
 end
     
 for i=1:length(algorithms)
-    algorithms(i).function(conn,CONSTANTS,algorithms(i).params,t);
+    featureTable=algorithms(i).function(conn,CONSTANTS,algorithms(i).params,t);
+    CellFeatures.writeFeatureTable(conn,featureTable);
 end
     

matlab/+CellFeatures/writeFeatureTable.m

+% should be cellID and 1 feature. write to db.
+function writeFeatureTable(conn,featureTable)
+
+if size(featureTable,2)~=2 || ~strcmp(featureTable.Properties.VariableNames{1},'cellID')
+    fprintf(2,'writeFeatureTable - invalid table format\n');
+    return;
+end
+    
+%     if isnan(ratio)
+%         cmd=['delete from tblCellFeatures where cellID=' num2str(cellIDs(i))];
+%         exec(conn,cmd);
+%     else
+%         ratio=round(ratio,4);
+%         cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
+%             num2str(cellIDs(i)) ')'];
+%         exec(conn,cmd);
+%         cmd=['update tblCellFeatures set ' params.name '=' num2str(ratio) ...
+%             ' where cellID=' num2str(cellIDs(i))];    
+%         exec(conn,cmd);
+%     end
+cmdDelete='delete from tblCellFeatures where cellID=?';
+psDelete=conn.handle.prepareStatement(cmdDelete);
+cmdInsert='insert or ignore into tblCellFeatures(cellID) values (?)';
+psInsert=conn.handle.prepareStatement(cmdInsert);
+cmdUpdate=['update tblCellFeatures set ' featureTable.Properties.VariableNames{2}...
+    '=? where cellID=?'];
+psUpdate=conn.handle.prepareStatement(cmdUpdate);
+
+for i=1:size(featureTable,1)
+    if isnan(featureTable{i,2})
+        psDelete.setDouble(1,featureTable{i,1});
+        Write.tryUpdateDB(psDelete);
+    else
+        psInsert.setDouble(1,featureTable{i,1});
+        Write.tryUpdateDB(psInsert);
+        % set feature value
+        psUpdate.setDouble(1,featureTable{i,2});
+        % set cellID
+        psUpdate.setDouble(2,featureTable{i,1});
+        Write.tryUpdateDB(psUpdate);
+    end
+end

matlab/+Ensemble/getNestedCells.m

@@ -61,6 +61,9 @@ for i=idx+1:length(ensembleCells)
                 if bucketMap(idxK)>0
                     idxK=bucketMap(idxK);
                 end
+                if idxK==idxDest
+                    continue
+                end
                 nestedCells{idxDest}=[nestedCells{idxDest},...
                     nestedCells{idxK}];
                 nestedCells{idxK}=[];

matlab/+Ensemble/setCellBackgroundEfficiency.m

-
-% sets the cells boundary efficiency feature. this is the one minus the average
-% along the a neighborhood of the boundary (or the vertices) of the percent 
-% of the gray threshold that the cell achieves. It is a measure of how dark
-% it is around the cell boundary. 1 is good (very dark boundary) 0 is bad
-% -- boundary very bright, indicative of cutting a cell in two
-function Cells=setCellBackgroundEfficiency(Cells,im,segParams)
-
-global DRAW
-
-% do a neighborhood min operator
-if is3D(im)
-    se=strel('sphere',3);
-else
-    se=strel('disk',3);
-end
-imBG=imdilate(im,se);
-
-for i=1:length(Cells)
-    if segParams.isPhase
-        bgEfficiency=1;
-    else
-        if is3D(im)
-            idx=sub2ind(size(im),Cells(i).pts(:,2),Cells(i).pts(:,1),Cells(i).pts(:,3));
-        else
-            idx=sub2ind(size(im),Cells(i).pts(:,2),Cells(i).pts(:,1));
-        end
-        [~,bgEfficiency]=graythresh(imBG(idx));
-    end
-    % ACK - efficiency or deficiency?
-    Cells(i).ensemble.backgroundEfficiency=1-bgEfficiency;
-end
-    

matlab/+Ensemble/setCellBoundaryEfficiency.m → matlab/+Ensemble/setCellEfficiency.m

@@ -4,7 +4,7 @@
 % of the gray threshold that the cell achieves. It is a measure of how dark
 % it is around the cell boundary. 1 is good (very dark boundary) 0 is bad
 % -- boundary very bright, indicative of cutting a cell in two
-function Cells=setCellBoundaryEfficiency(Cells,im,segParams)
+function Cells=setCellEfficiency(Cells,im,segParams)
 
 global DRAW
 
@@ -20,13 +20,10 @@ if isfield(segParams,'sensitivity')
 else
     sensitivity=0.5;
 end
-if sensitivity>1
-    lm=multithresh(im,sensitivity);
-    T=lm(end);
-else
+
 T=adaptthresh(im,sensitivity,'NeighborhoodSize',4*floor(size(im)/16)+1,'statistic','gaussian');
 % T=imdilate(T,se);
-end
+
 bw=logical(0*im);
 for i=1:length(Cells)
     idx=sub2ind(size(im),Cells(i).pts(:,2),Cells(i).pts(:,1));
@@ -39,22 +36,23 @@ for i=1:length(Cells)
         verts=max(verts,1);
         verts=min(verts,size(im));       
         idx=sub2ind(size(im),verts(:,1),verts(:,2),verts(:,3));
+        idxPts=sub2ind(size(im),Cells(i).pts(:,2),Cells(i).pts(:,1),Cells(i).pts(:,3));           
     else
         idx=sub2ind(size(im),Cells(i).surface(:,2),Cells(i).surface(:,1));
+        idxPts=sub2ind(size(im),Cells(i).pts(:,2),Cells(i).pts(:,1));
     end
     idx(isnan(idx))=[];    
-    imBGmu=median(imBG(idx));
-    if sensitivity>1
-        minFG=T;
-        imBGeff=1-imBGmu/minFG; %
+    % high number means boundary pixels are close to background
+    imBoundaryEff=1-mean(imBG(idx)-T(idx));
+    if segParams.isPhase
+        imBackgroundEff=1;
     else
-        imBGeff=1-median(imBG(idx)-T(idx));
+        % high number means cell pixels are brighter than background
+        imBackgroundEff=mean(im(idxPts)-T(idxPts));
     end
-    % clamp??? ac 10 14 2020 ??? ACK
-    %     imBGeff=min(imBGeff,1);
-    %     imBGeff=max(imBGeff,0);
     
     % how dark on average the surface neighborhood is (1 is good)
-    Cells(i).ensemble.boundaryEfficiency=imBGeff;
+    Cells(i).ensemble.boundaryEfficiency=imBoundaryEff;
+    Cells(i).ensemble.backgroundEfficiency=imBackgroundEff;
 end
     

matlab/+Segment/frameSegment.m

@@ -35,7 +35,7 @@ if 3==length(segParams.minimumRadius_um)
     return;
 end
 bEnsemble=isfield(args,'ensemble_minimumRadius_um');
-
+bEnsemble=true;
 global USE_CUDA
 USE_CUDA=segParams.useCuda&&~ismac && HIP.Cuda.DeviceCount>0;
 tic;
@@ -155,8 +155,7 @@ end
 
 % % on ensembleSeg
 if bEnsemble
-    Cells=Ensemble.setCellBoundaryEfficiency(Cells,im,segParams);
-    Cells=Ensemble.setCellBackgroundEfficiency(Cells,im,segParams);
+    Cells=Ensemble.setCellEfficiency(Cells,im,segParams);
     if DRAW
         for i=1:length(Cells)
             szBG=[ num2str(Cells(i).ensemble.backgroundEfficiency,1) ',' ...

matlab/+Segment/thresholdImages.m

@@ -107,6 +107,7 @@ else
     pixelIdxList={rp.PixelIdxList};
 end
 idxReset=find(convDef<1.1);
+if ~isempty(idxReset)
     pixelIdxReset=vertcat(pixelIdxList{idxReset});
     bwLoG(pixelIdxReset)=0;
-
+end

matlab/+Smooth/+Classify/goClassifyMitosis.m

@@ -2,6 +2,11 @@
 % pass in [conn,CONSTANTS] (from openDB) or strDB (full path to .LEVER)
 function classifyMitosis(conn,CONSTANTS)
 
+reseg=Read.getAlgorithms(conn,'reseg');
+if isempty(reseg) || ~strcmp(reseg.params.mitosisDetection,'true')
+    return;
+end
+
 mParams=Read.getMitosisParams(conn);
 targetTracks=Smooth.getTargetTracks(conn,5);
 strDB=conn.DataSource;