...
 
Commits (64)
......@@ -386,7 +386,7 @@ function getMCRpath()
return mcrPath;
// unix 2017b installed
mcrPath='/usr/local/MATLAB/R2017b/';
mcrPath='/usr/local/MATLAB/R2019b/';
if (fs.existsSync(mcrPath))
return mcrPath;
// unix 2017b MCR
......
......@@ -22,52 +22,6 @@ function anyFileInQ(processQ,leverFileList)
return null;
} // anyFileInQ
function batchResegAll(leverFileList)
{
const process=require('./process.js');
for (var i=0;i<leverFileList.length;i++) {
var leverCommand={leverFile:leverFileList[i],time:1,command:'resegAll',params:[1]};
process.addToQ(leverCommand);
}
} // batchResegAll
var batchTimer=null;
function monitorBatchList(leverFileList)
{
const process=require('./process.js');
if (null!==batchTimer) {
clearTimeout(batchTimer);
batchTimer=null;
}
if (undefined==leverFileList || 0==leverFileList.length)
return; // all done
for (var i=0;i<batchFilenames.length;i++) {
process.getLeverFileStatus(batchFilenames[i],function(leverFile,uiStatus){
if ('batchTrackPending'==uiStatus) {
var idx=batchFilenames.indexOf(leverFile);
if (idx>=0)
batchFilenames.splice(idx,1);
if (0==batchFilenames.length){
// reseg all
batchResegAll(leverFileList);
}
}
});
}
batchTimer=setTimeout(monitorBatchList,1000,leverFileList);
} // monitorBatchList
var batchFilenames=null;
function startBatchMonitor(leverFileList)
{
batchFilenames=[];
for (var i=0;i<leverFileList.length;i++)
batchFilenames.push(leverFileList[i]);
monitorBatchList(leverFileList);
} // startBatchMonitor
function batchSeg(leverFileList,fnCallback)
{
const process=require('./process.js');
......@@ -80,16 +34,12 @@ function batchSeg(leverFileList,fnCallback)
inUseFile=path.basename(inUseFile);
fnCallback(409,"leverFile :: "+inUseFile+" :: is already in the process Q -- cannot batch segment");
return;
}
let nComplete=0;
}
for (let i=0;i<leverFileList.length;i++) {
var leverDB=new LeverDB(leverFileList[i]);
leverDB.updateStatus('batch seg pending',function(){
leverDB.updateStatus('batchSeg pending',function(){
var leverCommand={leverFile:leverFileList[i],time:1,command:'batchSeg',params:[]};
process.addToQ(leverCommand);
nComplete++;
if (nComplete>=leverFileList.length)
startBatchMonitor(leverFileList);
});
}
fnCallback(200,"OK");
......
......@@ -11,7 +11,7 @@ function optimize(dbFile,bQuiet=false)
var cmd;
cmd='PRAGMA optimize';
lDB.exec(cmd,function(){
cmd='PRAGMA wal_checkpoint(full)';
cmd='PRAGMA journal_mode=WAL;PRAGMA wal_checkpoint(full)';
lDB.exec(cmd,function(){
lDB.close(function (error) {
if (null!=error)
......
......@@ -3,7 +3,7 @@ const fs = require('fs');
// use sha1 so it will match our git commit tags
const hash = crypto.createHash('sha1');
//const input = fs.createReadStream('..\\matlab\\+Segment\\FrameSegment.m');
//const input = fs.createReadStream('..\\matlab\\+Segment\\frameSegment.m');
if (3!=process.argv.length) {
console.error('USAGE: node goHash.js filename');
......
......@@ -1297,8 +1297,16 @@ class LeverDB {
fnCallback(vertices);
});
}
// set write ahead mode on the database -- makes concurrent access possible
setWAL(fnCallback) {
var sqlCmd = "PRAGMA journal_mode = WAL;pragma journal mode";
this._lDB.all(sqlCmd,function(err,rows){
fnCallback();
});
}
} // class leverDB
}
// read all the leverfiles in a folder, send one massive response
// unused, since it takes
// note - currently unused...
......
......@@ -701,6 +701,8 @@ function getLineageMap(rowsFamilies,rowsTracks,rowsPhenotypes)
// the phenotypes are indexed into constants.phenotypes (TBD). phenotype 0 is death.
for (i=0;i<rowsPhenotypes.length;i++) {
node=gLineageMap.get(rowsPhenotypes[i].trackID);
if (undefined==node)
continue;
if (undefined==node.phenotypes)
node.phenotypes=[];
if (undefined==node.phenotypes[rowsPhenotypes[i].time])
......
......@@ -64,7 +64,7 @@ function checkToProcessQ()
internalAddToQ(leverCommand[0]);
}
}
toProcessQtimer=setTimeout(checkToProcessQ,1000);
toProcessQtimer=setTimeout(checkToProcessQ,2000);
} // checkToProcessQ
var toProcessQtimer;
......@@ -112,16 +112,11 @@ function internalAddToQ(leverCommand)
// note - when batchSeg is active, we disallow subsequent batch segs
if (null!=batchSegFile && 'batchSeg'==leverCommand.command) {
addToProcessQ(leverCommand);
getLeverFileStatus(batchSegFile,function(leverFile,uiStatus){
if ('batchTrackPending'==uiStatus || 'idle'==uiStatus) {
batchSegFile=null;
}
});
return;
}
if ('batchSeg'==leverCommand.command) {
batchSegFile=leverCommand.leverFile;
}
}
// first, see if the leverFile is in the processQ
for (var i=0;i<inProcessQ.length;i++) {
if (undefined===inProcessQ[i].leverFile || null===inProcessQ[i].leverFile)
......@@ -232,6 +227,7 @@ function getUserID(leverCommand,fnCallback)
if (undefined!==fnCallback)
fnCallback(userID);
}
leverDB.close();
});
});
}
......@@ -241,6 +237,7 @@ function getUserID(leverCommand,fnCallback)
userID=rows[0].userID;
if (undefined!==fnCallback)
fnCallback(userID);
leverDB.close();
}
});
......@@ -256,7 +253,6 @@ function addToQ(leverCommand)
if (gEditID>Number.MAX_SAFE_INTEGER)
gEditID=1;
}
getUserID(leverCommand,function(userID){
if (undefined===leverCommand.params) {
leverCommand.params={};
......@@ -322,8 +318,13 @@ function addFiletoBatchSeg(leverCommand,nSlot)
});
inProcessQ[nSlot].childProc.stdout.on('data', (data) => {
// ack -- stdout?
ljsLog.log(' nSlot:'+nSlot+' :: '+data,1);
ljsLog.log(' nSlot:'+nSlot+' :: '+data,1);
var strout=data.toString();
var rexp=strout.match(/matlabPollDB pid=(?<pid>\d+)/);
if (null!=rexp) {
inProcessQ[nSlot].matlabPID=parseInt(rexp.groups.pid);
}
});
inProcessQ[nSlot].childProc.stderr.on('data', (data) => {
......@@ -331,11 +332,21 @@ function addFiletoBatchSeg(leverCommand,nSlot)
});
inProcessQ[nSlot].childProc.on('close', (code) => {
ljsLog.log('inProcessQ slot '+nSlot+' : exited with code ' + code);
var lf=inProcessQ[nSlot].leverFile;
inProcessQ[nSlot].childProc=null;
inProcessQ[nSlot].leverFile=null;
inProcessQ[nSlot].status="";
checkToProcessQ();
inProcessQ[nSlot].matlabPID=0; // done with the batch segment process
batchSegFile=null; // ok to start a new batch segment
// start a reseg on this file
// note - lf will be null if this process was just killed by the user. if so, don't try to keep going...
if (null!==lf) {
var leverCommand={leverFile:lf,time:1,command:'resegAll',params:[1]};
addToQ(leverCommand);
}
});
} // addFiletoBatchSeg
......@@ -344,7 +355,13 @@ function addFiletoBatchSeg(leverCommand,nSlot)
function addFileToProcess(leverCommand,nSlot)
{
if ("batchSeg"==leverCommand.command){
addFiletoBatchSeg(leverCommand,nSlot);
const LeverDB=require('./lever.js');
var leverDB=new LeverDB(leverCommand.leverFile);
leverDB.setWAL(function(){
leverDB.close(function(){
addFiletoBatchSeg(leverCommand,nSlot);
});
});
return;
}
const path=require('path');
......@@ -500,6 +517,7 @@ function removeID(editID)
function pauseOrKillAll(bKill,leverFile)
{
const spawn = require('child_process').spawn;
toProcessQ=[];
for (var i=0;i<inProcessQ.length;i++) {
......@@ -511,11 +529,21 @@ function pauseOrKillAll(bKill,leverFile)
if (bKill) {
gCommandMap.delete(inProcessQ[i].ID);
if (null!==inProcessQ[i].childProc)
inProcessQ[i].childProc.kill();
if (null!==inProcessQ[i].childProc) {
if (process.platform != "win32" && undefined!=inProcessQ[i].matlabPID) {
// matlab was launched on unix via a script
// it told us its pid right after launch
// we use that to kill the main matlab process here
// killing the script leaves the process running
process.kill(inProcessQ[i].matlabPID);
} else {
inProcessQ[i].childProc.kill();
}
}
inProcessQ[i].childProc=null;
inProcessQ[i].leverFile=null;
batchSegFile=null;
inProcessQ[nSlot].matlabPID=0;
}
else {
// pause gets put right into Q
......
......@@ -950,12 +950,12 @@ function drawHulls(tHulls)
}
// color by feature?
if (nFeature>=0) {
var fmax=gFeatureExtrema[nFeature][1];
var fmin=gFeatureExtrema[nFeature][0];
// scale so colors go from [0.2,1]
// puts it on [0,1]
var fval=1 - (tHulls[i][CONSTANTS.cellFeatures[0].name] - fmin)/((fmax-fmin));
// now put it on [0.1,1]
// var fmax=gFeatureExtrema[nFeature][1];
// var fmin=gFeatureExtrema[nFeature][0];
// // scale so colors go from [0.2,1]
// // puts it on [0,1]
var fval=tHulls[i][CONSTANTS.cellFeatures[0].name];
// // now put it on [0.1,1]
fval=0.9*fval+0.1;
color = getHTMLcolor([fval,fval,fval]);
}
......
......@@ -48,7 +48,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1
sp.NLM(4)=commandList{labindex+idxCommandList,'hpf'};
sp.alpha=commandList{labindex+idxCommandList,'alpha'};
sp.ensemble_minimumRadius_um=true;
nc=Segment.FrameSegment_texture('',1,CONSTANTS,sp);
nc=Segment.frameSegment('',1,CONSTANTS,sp);
else
nc=[];
end
......
......@@ -27,6 +27,15 @@ else
end
segParams.draw=false;
commandList=[];
if 1==segParams.wellRadius
% need to set the median image before we go into spmd, otherwise out of
% memory...
medianMask=Segment.getMedianMask(conn,CONSTANTS,segParams);
else
medianMask=[];
end
% make the command list
for t=1:CONSTANTS.imageData.NumberOfFrames
nc=[repmat(t,length(rgRadius),1),rgRadius'];
......@@ -68,7 +77,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1
sp=segParams;
sp.minimumRadius_um=commandList(labindex+idxCommandList,2);
if isempty(intersect(t,tSeg)) && t>0 && t<=CONSTANTS.imageData.NumberOfFrames
nc=Segment.FrameSegment_texture('',t,CONSTANTS,sp);
nc=Segment.frameSegment('',t,CONSTANTS,sp,medianMask);
else
if ~isempty(intersect(t,tSeg))
fprintf(1,' found existing seg, t=%d, idxCommandList=%d\n',t,idxCommandList+labindex);
......
......@@ -19,11 +19,16 @@ for ff=1:length(flist)
Batch.batchSegment(strDB,[],false);
end
if ispc
nodecmd='node.exe ';
else
nodecmd='node ';
end
ljsPath=getLeverjsPath('resegAndTrackDataset.js');
if bTrack
parfor ff=1:length(flist)
strDB=fullfile(flist(ff).folder,flist(ff).name);
system(['node.exe ' ljsPath ' "' ...
system([nodecmd ' ' ljsPath ' "' ...
strDB '" resegAll 1 "{}" ' ]);
end
end
......
......@@ -5,7 +5,7 @@ exec(conn,['ALTER TABLE tblCellFeatures ADD COLUMN ' featureName ' NUMBER']);
nf.name=featureName;
nf.color=color;
if ~isfield(CONSTANTS,'cellFeatures')
if ~isfield(CONSTANTS,'cellFeatures') || isempty(CONSTANTS.cellFeatures)
CONSTANTS.cellFeatures=nf;
else
if ~any(cell2mat(strfind({CONSTANTS.cellFeatures.name},nf.name)))
......
function membraneRatio(conn,CONSTANTS,args,t)
params.channel=1;
params.name='membraneRatio';
params.color="#ff0000";
% check for param override
params=Segment.getParams(params,args);
% segmentation params -- for min radius
segParams=Read.getSegmentationParams(conn);
bw=Segment.getMembrane(CONSTANTS,params.channel,t);
if isempty(bw) || any(isnan(bw(:))) || all(bw(:)==0)
return;
end
if is3D(bw)
fprintf(2,'membraneRatio untested in 3d - not computing\n');
return;
end
Cells=Read.getCellsTime(conn,t);
for i=1:length(Cells)
idxSurface=sub2ind(size(bw),Cells(i).surface(:,2),Cells(i).surface(:,1));
idxSurface(isnan(idxSurface))=[];
ratio = length(find(bw(idxSurface)))/length(idxSurface);
cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
num2str(Cells(i).cellID) ')'];
exec(conn,cmd);
cmd=['update tblCellFeatures set ' params.name '=' num2str(ratio) ...
' where cellID=' num2str(Cells(i).cellID)];
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
end
ratio=median(imOuter)/median(imInner);
function [Cells,nestedCells]=ensembleSegment(conn, t,CONSTANTS,segParams)
function [Cells,nestedCells]=ensembleSegment(conn, t,CONSTANTS,segParams,medianMask)
ensembleDraw=false;
if segParams.draw==true
......@@ -37,7 +37,7 @@ parfor i=1:length(rgRadius)
% note as of matlab2019b, we pass [] for the conn (1st) argument as the
% param is unused as long as CONSTANTS are provided, and even passing
% it around in a parfor generates a warning
cx=Segment.FrameSegment_texture([],t,CONSTANTS,args);
cx=Segment.frameSegment([],t,CONSTANTS,args,medianMask);
ensembleCells{i}=cx;
end
segParams.draw=ensembleDraw;
......
......@@ -9,12 +9,12 @@ nctime=tic();
% initialize the bucket with the cells from the maximal radius
% segmentation
idx=1;
while isempty(ensembleCells{idx}) && idx<length(ensembleCells)
while isempty(ensembleCells{idx}) && idx<=length(ensembleCells)
idx=idx+1;
end
if idx==length(ensembleCells)
if idx>length(ensembleCells)
% only zero or one non-empty bucket -- nothing to nest
nestedCells=ensembleCells{idx}; % could be [], or just cells from smallest radius
nestedCells={};
return;
end
for i=1:length(ensembleCells{idx})
......@@ -30,6 +30,9 @@ end
% the two buckets
for i=idx+1:length(ensembleCells)
iCells=ensembleCells{i};
if isempty(iCells)
continue
end
p=gcp('nocreate');
if isempty(p)
[iBucketDest,iCells]=getBucketDest(iCells,rgRadius,nestedCells,szIm,i);
......
......@@ -15,18 +15,26 @@ else
se=strel('disk',3);
end
imBG=imerode(im,se);
if isfield(segParams,'sensitivity')
sensitivity=segParams.sensitivity;
else
sensitivity=0.5;
end
T=adaptthresh(im,sensitivity,'NeighborhoodSize',4*floor(size(im)/16)+1,'statistic','gaussian');
T=imdilate(T,se);
for i=1:length(Cells)
if is3D(im)
verts=round(Cells(i).verts);
verts=[verts(:,2),verts(:,1),verts(:,3)];
verts=max(verts,1);
verts=min(verts,size(im));
idx=sub2ind(size(im),verts(:,2),verts(:,1),verts(:,3));
idx=sub2ind(size(im),verts(:,1),verts(:,2),verts(:,3));
else
idx=sub2ind(size(im),Cells(i).surface(:,2),Cells(i).surface(:,1));
end
imBGmu=mean(imBG(idx));
minFG=graythresh(im)*segParams.alpha(1);
idx(isnan(idx))=[];
imBGmu=median(imBG(idx));
minFG=max(T(idx));
imBGeff=1-imBGmu/minFG; %
imBGeff=min(imBGeff,1);
......
% called by Segmente.FrameSegment_texture to set the
% called by Segmente.frameSegment to set the
% area, convexArea, circularity on each cell
function newCell=setCellFeatures(newCell,RP,n,segParams)
......
......@@ -55,7 +55,7 @@ for t=1:CONSTANTS.imageData.NumberOfFrames
if length(idx)<5
continue
end
newCell=Segment.FrameSegment_3D_create(idx,size(bwMeta),metaChannel,t);
newCell=Segment.frameSegment_create_3D(idx,size(bwMeta),metaChannel,t);
metaCells=[metaCells newCell];
end
end
......
......@@ -2,12 +2,11 @@
UISERVER='http://localhost:4444/'
movieFile='pos1_rotate_seg.mp4';
movieFile='9drotation.mp4';
v=VideoWriter(movieFile,'MPEG-4');
v.FrameRate=10;
open(v)
ui=webread([UISERVER 'ui']);
ui.sidebar='none';
ui.webToolbar='none';
......@@ -21,19 +20,6 @@ webwrite([UISERVER 'ui'],options,ui);
time=1;
webwrite([UISERVER 'time/' num2str(time)] ,'');
% wait for mr render to complete
bComplete=false;
while ~bComplete
bComplete=webread([UISERVER 'drawComplete']);
if ~bComplete
pause(0.5)
end
end
im=webread([UISERVER 'screenCap']);
imagesc(im)
axis image
axis off
view=webread([UISERVER 'view']);
options = weboptions('MediaType','application/json');
theta=1/pi;
......@@ -42,11 +28,8 @@ theta=1/pi;
% rz=[cos(theta),-sin(theta),0,0;sin(theta),cos(theta),0,0;0,0,1,0;0,0,0,1];
for theta=0:0.1:2*pi
ry=[cos(theta), 0, sin(theta),0;0,1,0,0;-sin(theta),0,cos(theta),0;0,0,0,1];
%worldRot=reshape(view.worldRot,4,4);
worldRot=ry;
% worldRot=worldRot';
view.worldRot=worldRot(:);
% view.zoom=0.5;
webwrite([UISERVER 'view'],options,view);
bComplete=false;
while ~bComplete
......@@ -56,9 +39,9 @@ for theta=0:0.1:2*pi
end
end
im=webread([UISERVER 'screenCap']);
% imagesc(im)
% axis image
writeVideo(v,im);
imagesc(im)
axis image
% writeVideo(v,im);
drawnow
fprintf(1,'theta=%f complete\n',theta);
end
......
......@@ -25,7 +25,7 @@ end
srcTracks=cell2mat(q);
% tracks that start in this frame
cmd=['select cellID,time from tblCells where cellID=trackID and time>='...
num2str(tReseg-5) ' and time<=' num2str(tReseg+1) ...
num2str(tReseg-2) ' and time<=' num2str(tReseg+1) ...
' and cellID not in (select cellID_dst from tblManualEdits)'];
q=ljsFetch(conn,cmd);
if isempty(q)
......
......@@ -15,7 +15,7 @@ editQ=ljsFetch(conn,cmd);
if ~isempty(editQ)
idxUndo = cellfun(@(x) strcmp(x,'undo'),editQ(:,5));
idxUndone = [editQ{idxUndo,1}]'
idxUndone = [editQ{idxUndo,1}]';
for i=1:size(editQ,1)
if idxUndo(i)
continue;
......
......@@ -84,7 +84,7 @@ end
segParams=Read.getSegmentationParams(conn);
ptsReplicate_xy=Segment.FrameSegment_texture_getSplitPixels(conn,CONSTANTS,...
ptsReplicate_xy=Segment.frameSegment_getSplitPixels(conn,CONSTANTS,...
idxPixels, channel, time, segParams);
try
......@@ -106,7 +106,7 @@ for kpts=1:max(idx)
return;
end
ptsIdx=find(idx==kpts);
newCell=Segment.FrameSegment_texture_create(idxPixels(ptsIdx),size(bw),...
newCell=Segment.frameSegment_create(idxPixels(ptsIdx),size(bw),...
channel,time,conn,segParams);
newCell.segCC=Q{6};
Cells=[Cells newCell];
......
function [Cells,nestedCells]=FrameSegment_texture(conn, t,CONSTANTS,args)
Cells=[];
nestedCells=[];
if nargin<3
CONSTANTS=Read.getConstants(conn);
end
segParams=Segment.getDefaultSegParams();
if nargin==4
segParams=Segment.getParams(segParams,args);
end
global DRAW
DRAW=segParams.draw;
% ensemble segmentations
if 3==length(segParams.minimumRadius_um)
segParams.ensemble_minimumRadius_um=segParams.minimumRadius_um;
[Cells,nestedCells]=Ensemble.ensembleSegment(conn, t,CONSTANTS,segParams);
return;
end
bEnsemble=isfield(args,'ensemble_minimumRadius_um');
bEnsemble=true;
global USE_CUDA
USE_CUDA=segParams.useCuda&&~ismac && HIP.Cuda.DeviceCount>0;
tic;
resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
min_radius_pixels = segParams.minimumRadius_um ./ resolution_um;
chan=segParams.channels(1); % default channel
[im,qLog,imLog]=Segment.getImages(CONSTANTS,t,chan,segParams);
if all(im(:)==0)
ljsLog(sprintf('frameSegment_texture: empty image at time=%d -- skipping\n',t));
return;
end
for idChannel=2:length(segParams.channels)
[~,channel_qLog]=Segment.getImages(CONSTANTS,t,segParams.channels(idChannel),segParams);
qLog=qLog+channel_qLog;
end
if is3D(im)
min_area_pixels = 4/3*pi* mean(min_radius_pixels)^3;
else
min_radius_pixels=min_radius_pixels(1);
min_area_pixels = min_radius_pixels^2 * pi;
end
min_area_pixels = max(1,round(min_area_pixels));
if is3D(im)
if USE_CUDA
entropyRadius=round(min_radius_pixels);
kernel=HIP.MakeEllipsoidMask(entropyRadius);
imTexture=HIP.EntropyFilter(im,kernel,[]);
else
imTexture=im;
end
else
% 2D
if USE_CUDA
entropyRadius=round(min_radius_pixels);
kernel=HIP.MakeEllipsoidMask([entropyRadius,entropyRadius,0]);
imTexture=HIP.EntropyFilter(im,kernel,[]);
else
se=strel('square',2*floor(min_radius_pixels/2)+1);
imTexture=entropyfilt(double(im),getnhood(se));
end
imTexture=mat2gray(imTexture);
end
if DRAW
imRaw=MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',segParams.channels(1),'timeRange',[t t], 'outType','single','prompt',false);
imc=mat2gray(imRaw);
imq=mat2gray(imTexture);
if is3D(imc)
imc=max(imc,[],3);
imq=max(imq,[],3);
end
% hold off;imagesc([imc imq]);hold on
hold off;imagesc([imc]);colormap(gray);brighten(0.5);hold on
drawnow
end
[bw,bwIntensity,nLevels]=Segment.thresholdImages(im,imTexture,segParams,...
min_radius_pixels,min_area_pixels);
sz=sprintf('frameSegmentTexture: min_radius_pixels=%s nLevels detected at %d\n',mat2str(min_radius_pixels,3),nLevels);
ljsLog(sz,3);
bw=bw|imfill(bw,'holes');
if segParams.wellRadius>0
% mask all points at r>radius
center=size(im)/2;
rMask=ones(size(im));
[r,c]=find(rMask);
radius= sqrt( (r-center(1)).^2+(c-center(2)).^2 );
rMask=reshape(radius,size(im));
bwMask=logical(rMask<segParams.wellRadius);
bw(~bwMask)=0;
end
if ~segParams.bCytoplasmic
bw=removeOuterRadius(bw,CONSTANTS,min_radius_pixels);
end
qLog(~bw)=0;
% bw2 is all the pixels in our foreground image. after we segReduce, some
% of those pixels will be discarded in resolving the underlying cells.
% allocateShake (below) assigns each discarded pixel to it's closest
% touching segmentation using a morphological region fill
bw2=bwareaopen(bw,min_area_pixels);
% tic
bw=Segment.segReduce(bw,qLog,min_area_pixels,min_radius_pixels);
% tReduce=toc;
% fprintf(1,'segReduce time = %f\n',tReduce);
bw=bw|imfill(bw,'holes');
if DRAW
if is3D(bw)
bounds=bwboundaries(max(bw,[],3));
else
bounds=bwboundaries(bw);
end
for i=1:length(bounds)
plot(bounds{i}(:,2),bounds{i}(:,1),'-r')
end
end
if isempty(bw)
bw=bw2;
end
preShakeL=bwlabeln(bw);
% ok. bw is the texture image, with connected components split up by
% sub-regions via segReduce (using quantized LOG image qLog).
% bwIntensity is the otsu thresholded intensity image. if a 2 element
% segParams.alpha is provided, the second element is the weight for the
% intensity otsu. otherwise, the single first element of segParams.alpha
% weights both texture and intensity thresholding.
% now, here we prepare to 'allocateShake' -- assigning texture foreground
% pixels removed in segReduce back to the nearest (topographically by qLog)
% cell region from bw.
% new in August 2019 we don't assign back pixels that were not foreground
% in either bw (our cell image from segreduce) or bwIntensity (our
% intensity thresholded image)
bw2 = bw2 & (bw | bwIntensity);
% tic
[L,num]=Segment.allocateShake(bw,bw2,qLog,min_radius_pixels);
% tShake=toc;
% fprintf(1,'allocateShake: t=%f\n',tShake);
% [L,num]=bwlabeln(bw);
[L2]=bwlabeln(bw2); % used to mark original cc for segCC field
preShakeL=L;
preShakeL(~bw)=0;
if bEnsemble
if size(bw,3)>1
rp=regionprops3(L,'Volume','ConvexVolume');
else
rp=regionprops(L,'Area','ConvexArea','Eccentricity');
end
else
rp=[];
end
% final CC pass
for n=1:num
idx=find(L==n);
if length(idx)<min_area_pixels
continue
end
newCell=Segment.FrameSegment_texture_create(idx,size(bw),chan,t);
segL=unique(L2(idx));
segL(0==segL)=[];
if length(segL)>1
segL=segL(1);
end
newCell.segCC=segL;
newCell.idxPreShake=find(preShakeL==n);
if bEnsemble
newCell=Ensemble.setCellFeatures(newCell,rp,n,segParams);
end
if DRAW && ~is3D(im)
for i=1:length(newCell)
nc=newCell(i);
cmap=hsv(255);
cidx=1+round(rand()*254);
plot(nc.surface(:,1),nc.surface(:,2),'color',cmap(cidx,:),'linewidth',2);
% drawnow
end
end
Cells=[Cells newCell];
end
% set features for parasite test
for i=1:length(Cells)
if is3D(bwIntensity)
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
Cells(i).pfg=length(find(bwIntensity(idx)))/size(Cells(i).pts,1);
% less then 5% foreground are parasites
if DRAW && Cells(i).pfg<0.05
plot(Cells(i).surface(:,1),Cells(i).surface(:,2),'color','k','linewidth',3);
end
end
if ~isempty(Cells)
% if 5% or less of pixels are foreground, delete
idxParasite=find( [Cells.pfg]<0.05 );
Cells(idxParasite)=[];
end
% on ensembleSeg
if bEnsemble
Cells=Ensemble.setCellBoundaryEfficiency(Cells,im,segParams);
Cells=Ensemble.setCellBackgroundEfficiency(Cells,im,segParams);
if DRAW
for i=1:length(Cells)
szBG=[ num2str(Cells(i).ensemble.backgroundEfficiency,1) ',' ...
num2str(Cells(i).ensemble.boundaryEfficiency,1)];
% text(Cells(i).centroid(1),Cells(i).centroid(2),szBG,...
% 'background','w','fontsize',8);
end
end
end
tElapsed=toc;
sz=sprintf('segmented frame %d found %d hulls, useCuda=%d, time=%f\n',t,length(Cells),USE_CUDA,tElapsed);
ljsLog(sz,3);
if DRAW
title([num2str(tElapsed,'%0.0f') 'seconds']);
drawnow
end
4;
function bw=removeOuterRadius(bw,CONSTANTS,min_radius_pixels)
global USE_CUDA
if is3D(bw)
% resize bw so that pixels are isometric -- um in each dimension
scale=CONSTANTS.imageData.PixelPhysicalSize./min(CONSTANTS.imageData.PixelPhysicalSize);
newSize=size(bw).*scale;
% each pixel in new image will have a physical size
% min(CONSTANTS.imageData.PixelPhysicalSize) cube on each side
bwr=imresize3(mat2gray(bw),newSize);
bwd=bwdist(~bwr);
bwd=imresize3(bwd,size(bw));
bw(bwd<max(round(min_radius_pixels)))=0;
% now a closing to remove any 'islands' between cells
if USE_CUDA
entropyRadius=round(min_radius_pixels);
kernel=HIP.MakeEllipsoidMask(entropyRadius);
bw=HIP.Opener(bw,kernel,1,[]);
end
else
se=strel('disk',round(min_radius_pixels(1)));
% remove the filter overhang
bw=imerode(bw,se);
% now remove the 'islands'
bw=imopen(bw,se);
end
function [Cells,nestedCells]=frameSegment(conn, t,CONSTANTS,args, medianMask)
Cells=[];
nestedCells=[];
if nargin<3
CONSTANTS=Read.getConstants(conn);
end
segParams=Segment.getDefaultSegParams();
if exist('args','var')
segParams=Segment.getParams(segParams,args);
end
if isdeployed()
segParams.draw=false;
end
global DRAW
DRAW=segParams.draw;
if (1==segParams.wellRadius) && ~exist('medianMask','var')
medianMask=Segment.getMedianMask(conn,CONSTANTS,segParams);
else
if ~exist('medianMask','var')
medianMask=[];
end
end
% ensemble segmentations
if 3==length(segParams.minimumRadius_um)
segParams.ensemble_minimumRadius_um=segParams.minimumRadius_um;
[Cells,nestedCells]=Ensemble.ensembleSegment(conn, t,CONSTANTS,segParams,medianMask);
return;
end
bEnsemble=isfield(args,'ensemble_minimumRadius_um');
bEnsemble=true;
global USE_CUDA
USE_CUDA=segParams.useCuda&&~ismac && HIP.Cuda.DeviceCount>0;
tic;
resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
min_radius_pixels = segParams.minimumRadius_um ./ resolution_um;
chan=segParams.channels(1); % default channel
[im,imLog]=Segment.getImages(CONSTANTS,t,chan,segParams,medianMask);
if all(im(:)==0)
ljsLog(sprintf('frameSegment: empty image at time=%d -- skipping\n',t));
return;
end
if is3D(im)
min_area_pixels = 4/3*pi* mean(min_radius_pixels)^3;
else
min_radius_pixels=min_radius_pixels(1);
min_area_pixels = min_radius_pixels^2 * pi;
end
min_area_pixels = max(1,round(min_area_pixels));
if DRAW
imRaw=MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',segParams.channels(1),'timeRange',[t t], 'outType','single','prompt',false);
imc=mat2gray(imRaw);
if is3D(imc)
imc=max(imc,[],3);
end
hold off;imagesc([imc]);colormap(gray);axis image;brighten(0.5);hold on
drawnow
end
[bw,bwLog]=Segment.thresholdImages(im,imLog,segParams, min_radius_pixels,min_area_pixels,...
medianMask);
if segParams.membraneChannel>0
bwMembrane=Segment.getMembrane(CONSTANTS,segParams.membraneChannel,t);
if ~isempty(bwMembrane)
bwLog=bwLog|bwMembrane;
end
end
% bw2 is all the pixels in our foreground image. after we mask, some
% of those pixels will be discarded in resolving the underlying cells.
% allocateShake (below) assigns each discarded pixel to it's closest
% touching segmentation using a morphological region fill
bw2=bw;
bw=Segment.segReduce(bw,bwLog,segParams,min_area_pixels,min_radius_pixels);
if ~segParams.isPhase
bw=bw|imfill(bw,'holes');
end
[L,num]=bwlabeln(bw);
if ~is3D(im)
assignRadius=2*min_radius_pixels(1);
else
assignRadius=min(min_radius_pixels).*2;
end
[d,dNN]=bwdist(bw);
bwAssign=bw2&~bw;
bwAssign(d>assignRadius)=0;
idxAssign=find(bwAssign);
L(idxAssign)=L(dNN(idxAssign));
num=max(L(:));
% tic
% [L,num]=Segment.allocateShake(bw,bw2,qLog,min_radius_pixels);
[L2]=bwlabeln(bw2); % used to mark original cc for segCC field
if bEnsemble
if size(bw,3)>1
rp=regionprops3(L,'Volume','ConvexVolume');
else
rp=regionprops(L,'Area','ConvexArea','Eccentricity');
end
else
rp=[];
end
d=bwdist(~logical(L));
% final CC pass
for n=1:num
idx=find(L==n);
if length(idx)<2*min_area_pixels || max(d(idx))<min_radius_pixels/2
continue
end
newCell=Segment.frameSegment_create(idx,size(bw),chan,t);
segL=unique(L2(idx));
segL(0==segL)=[];
if length(segL)>1
segL=segL(1);
end
newCell.segCC=segL;
if bEnsemble
newCell=Ensemble.setCellFeatures(newCell,rp,n,segParams);
end
if DRAW && ~is3D(im)
for i=1:length(newCell)
nc=newCell(i);
cmap=hsv(255);
cidx=1+round(rand()*254);
plot(nc.surface(:,1),nc.surface(:,2),'color',cmap(cidx,:),'linewidth',1);
% drawnow
end
end
Cells=[Cells newCell];
end
% % intensity threshold for parasites
% if segParams.isPhase>0
% % bright interior phase
% fgRatio=0.05; % intensity model needed to keep phase parasites down
% elseif segParams.isPhase<0
% % dark interior phase
% fgRatio=0.05; % less intensity driven
% else
% % not phase
% fgRatio=0.05; % less intensity driven
% end
% % set features for parasite test
% for i=1:length(Cells)
% if is3D(bw)
% 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
% Cells(i).pfg=length(find(bwIntensity(idx)))/size(Cells(i).pts,1);
% % less then 5% foreground are parasites
% if DRAW
% % text(Cells(i).centroid(1),Cells(i).centroid(2),num2str(Cells(i).pfg,2),'color','w');
% if Cells(i).pfg<fgRatio
% plot(Cells(i).surface(:,1),Cells(i).surface(:,2),'color','k','linewidth',3);
% end
% end
% end
%
% if ~isempty(Cells)
% % if 5% or less of pixels are foreground, delete
% idxParasite=find( [Cells.pfg]<fgRatio );
% Cells(idxParasite)=[];
% end
%
% % on ensembleSeg
if bEnsemble
Cells=Ensemble.setCellBoundaryEfficiency(Cells,im,segParams);
Cells=Ensemble.setCellBackgroundEfficiency(Cells,im,segParams);
if DRAW
for i=1:length(Cells)
szBG=[ num2str(Cells(i).ensemble.backgroundEfficiency,1) ',' ...
num2str(Cells(i).ensemble.boundaryEfficiency,1)];
% text(Cells(i).centroid(1),Cells(i).centroid(2),szBG,...
% 'background','w','fontsize',8);
end
end
end
tElapsed=toc;
sz=sprintf('segmented frame %d found %d hulls, useCuda=%d, time=%f\n',t,length(Cells),USE_CUDA,tElapsed);
ljsLog(sz,3);
if DRAW
title([num2str(tElapsed,'%0.0f') 'seconds']);
drawnow
end
4;
function bw=removeOuterRadius(bw,CONSTANTS,min_radius_pixels)
global USE_CUDA
if is3D(bw)
% resize bw so that pixels are isometric -- um in each dimension
scale=CONSTANTS.imageData.PixelPhysicalSize./min(CONSTANTS.imageData.PixelPhysicalSize);
newSize=size(bw).*scale;
% each pixel in new image will have a physical size
% min(CONSTANTS.imageData.PixelPhysicalSize) cube on each side
bwr=imresize3(mat2gray(bw),newSize);
bwd=bwdist(~bwr);
bwd=imresize3(bwd,size(bw));
bw(bwd<max(round(min_radius_pixels)))=0;
% now a closing to remove any 'islands' between cells
if USE_CUDA
entropyRadius=round(min_radius_pixels);
kernel=HIP.MakeEllipsoidMask(entropyRadius);
bw=HIP.Opener(bw,kernel,1,[]);
end
else
se=strel('disk',round(min_radius_pixels(1)));
% remove the filter overhang
bw=imerode(bw,se);
% now remove the 'islands'
bw=imopen(bw,se);
end
function newCell=FrameSegment_texture_create(idxPixels, imSize, channel, time, ...
function newCell=frameSegment_create(idxPixels, imSize, channel, time, ...
conn, params)
if 3==length(imSize) && imSize(3)>1
% 3D
newCell=Segment.FrameSegment_3D_create(idxPixels, imSize, channel, time);
newCell=Segment.frameSegment_create_3D(idxPixels, imSize, channel, time);
return;
end
......
function newCell=FrameSegment_3D_create(idxPixels, imSize, channel, time, ...
conn, params)
bwn=zeros(imSize);
bwn(idxPixels)=1;
PAD=5;
% padd with extra zeros around the outside so our mesh covers boundaries
bpad = zeros(size(bwn)+10);
bpad(PAD+1:end-PAD,PAD+1:end-PAD,PAD+1:end-PAD)=bwn;
bwn = bpad;
% scale=[1,1,1];
scale=max(imSize)./imSize/8;
scale=max(min(scale,0.5),.01);
%
resize=round(scale.*size(bwn));
bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize));
if size(verts,1)<10
resize=round(2*scale.*size(bwn));
bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize));
end
norms = isonormals(bwResize,verts);
% rescale verts back to
unscale=(size(bwn)./size(bwResize));
% put unscale on x,y instead of r,c
unscale=[unscale(2),unscale(1),unscale(3)];
verts=(verts).*unscale-0.5*unscale;
edges=Segment.MakeEdges(faces);
% make edges and faces zero indexed
edges=edges-1;
faces=faces-1;
% subtract extract for padding
verts=verts-PAD;
maxRad = verts-repmat(mean(verts,1),size(verts,1),1);
maxRad=sum(abs(maxRad),2);
maxRad=max(maxRad);
newCell=[];
newCell.time=time;
newCell.centroid=[mean(verts,1)]-1;
newCell.edges=edges;
newCell.faces=faces;
newCell.verts=verts;
newCell.normals=norms;
xyzPts=[];
% note that xyzPts are on the correct range, as idxPixels is unpadded...
[xyzPts(:,2),xyzPts(:,1),xyzPts(:,3)]=ind2sub(imSize,idxPixels);
newCell.pts=uint16(xyzPts);
newCell.maxRadius=maxRad;
newCell.channel=channel;
function newCell=FrameSegment_create_3D(idxPixels, imSize, channel, time, ...
conn, params)
bwn=zeros(imSize);
bwn(idxPixels)=1;
PAD=5;
% padd with extra zeros around the outside so our mesh covers boundaries
bpad = zeros(size(bwn)+10);
bpad(PAD+1:end-PAD,PAD+1:end-PAD,PAD+1:end-PAD)=bwn;
bwn = bpad;
% scale=[1,1,1];
scale=max(imSize)./imSize/8;
scale=max(min(scale,0.5),.01);
%
resize=round(scale.*size(bwn));
bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize));
if size(verts,1)<10
resize=round(2*scale.*size(bwn));
bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize));
end
norms = isonormals(bwResize,verts);
% rescale verts back to
unscale=(size(bwn)./size(bwResize));
% put unscale on x,y instead of r,c
unscale=[unscale(2),unscale(1),unscale(3)];
verts=(verts).*unscale-0.5*unscale;
edges=Segment.MakeEdges(faces);
% make edges and faces zero indexed
edges=edges-1;
faces=faces-1;
% subtract extract for padding
verts=verts-PAD;
maxRad = verts-repmat(mean(verts,1),size(verts,1),1);
maxRad=sum(abs(maxRad),2);
maxRad=max(maxRad);
newCell=[];
newCell.time=time;
newCell.centroid=[mean(verts,1)]-1;
newCell.edges=edges;
newCell.faces=faces;
newCell.verts=verts;
newCell.normals=norms;
xyzPts=[];
% note that xyzPts are on the correct range, as idxPixels is unpadded...
[xyzPts(:,2),xyzPts(:,1),xyzPts(:,3)]=ind2sub(imSize,idxPixels);
newCell.pts=uint16(xyzPts);
newCell.maxRadius=maxRad;
newCell.channel=channel;
function ptsReplicate_xy=FrameSegment_texture_getSplitPixels(conn,CONSTANTS,idxPixels, channel, time, ...
function ptsReplicate_xy=frameSegment_getSplitPixels(conn,CONSTANTS,idxPixels, channel, time, ...
params)
segParams=Segment.getDefaultSegParams();
......
......@@ -2,19 +2,16 @@ function segParams=getDefaultSegParams()
segParams=[];
segParams.draw=false;
segParams.minimumRadius_um=[1.5,0.25,2.5];
segParams.channels=1; %[1CONSTANTS.imageData.NumberOfChannels];
segParams.channels=1;
segParams.membraneChannel=-1;
segParams.isPhase=false;
segParams.wellRadius=-1;
segParams.useCuda=true;
segParams.nCores=feature('numcores'); % for ensemble seg only
segParams.storeEnsemble=true;
segParams.alpha=1.0;
segParams.sensitivity=0.5;
segParams.denoise=true;
segParams.bCytoplasmic=false;
segParams.bCircular=false;
segParams.brightLevels=-1;
segParams.clonalMotionComp=false;
% advanced params
%
% NLM: h,search_window_radius,neighborhood_radius,high_pass_filter
segParams.NLM=[0.1,12,1,false];
function [im,qLog,imLogRaw]=getImages(CONSTANTS,t,channel,segParams)
function [im,imLog]=getImages(CONSTANTS,t,channel,segParams,medianMask)
if ~exist('bFilter','var')
bFilter=false;
......@@ -8,17 +8,22 @@ resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
min_radius_pixels = segParams.minimumRadius_um ./ resolution_um;
global USE_CUDA
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',channel,'timeRange',[t t], 'outType','single','prompt',false);
if 1==segParams.wellRadius && exist('medianMask','var') && ~isempty(medianMask)
im=im-medianMask;
im=mat2gray(im);
end
if isempty(im)
qLog=[];
imLogRaw=[];
imLog=[];
return