...
 
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 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
end
imRaw=im;
im=denoise(im,segParams);
if 1==nargout
return
end
if is3D(im)
% 3D
if 1==nargout
......@@ -26,43 +31,42 @@ if is3D(im)
end
if USE_CUDA
logRadius=min_radius_pixels.*(1/sqrt(3));
imLog=HIP.LoG(im,logRadius,[]);
imLog=HIP.LoG(imRaw,logRadius,[]);
else
% use difference of gaussian approximation
logRadius=min_radius_pixels;
imd1=imgaussfilt3(im,sqrt(2).*logRadius);
imd2=imgaussfilt3(im,logRadius./sqrt(2));
imd1=imgaussfilt3(imRaw,sqrt(2).*logRadius);
imd2=imgaussfilt3(imRaw,logRadius./sqrt(2));
imLog=imd1-imd2;
end
imLogRaw=imLog;
else
% 2D
if USE_CUDA
logRadius=(1/sqrt(2)).*min_radius_pixels;
logRadius(2)=logRadius(1);
logRadius(3)=0;
imLog=HIP.LoG(im,logRadius,[]);
imLog=HIP.LoG(imRaw,logRadius,[]);
else
if length(min_radius_pixels)>1
logRadius=0.5*min(min_radius_pixels(1:2));
else
logRadius=0.5*min_radius_pixels ;
end
szFilter=round(size(im)/3);
szFilter=round(size(imRaw)/3);
h=fspecial('log',szFilter,logRadius);
imLog=imfilter(im,h,'replicate');
imLog=imfilter(imRaw,h,'replicate');
end
end
imLogRaw=imLog;
imLog=mat2gray(imLog);
lm=multithresh(imLog,15);
qLog=imquantize(imLog,lm);
if segParams.isPhase
qLog=max(qLog(:))+1-qLog;
end
% note - DO NOT mat2gray imLog here. Later processing can use the negative
% or positive portions for inside/outside info. Once we mat2gray that's
% gone.
% LEFT HERE AS CAUTIONARY TALE imLog=mat2gray(imLog);
%
4;
function im=denoise(im,segParams)
if false==segParams.denoise
im=mat2gray(im);
return
......@@ -88,7 +92,7 @@ if is3D(im)
imx=imx-HIP.Cuda.Gaussian(im,szFilter,1,[]);
end
else
imx = medfilt3(im)-imgaussfilt3(im,szFilter);
imx = medfilt3(im);
end
else
% 2D
......@@ -98,10 +102,10 @@ else
imx=imx-HIP.Cuda.Gaussian(im,[szFilter,1],1,[]);
end
else
imx=medfilt2(im)-imgaussfilt(im,szFilter(1:2));
imx=imnlmfilt(im);
end
end
imx=max(imx,0);
im=mat2gray(imx);
......
function medianMask=getMedianMask(conn,CONSTANTS,segParams)
global medianMask
global datasetName
if isempty(datasetName) || isempty(medianMask) || ...
~strcmp(datasetName,fullfile(CONSTANTS.imageData.imageDir,CONSTANTS.imageData.DatasetName))
imMask = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',segParams.channels(1),'timeRange',[1 CONSTANTS.imageData.NumberOfFrames], 'outType','single','prompt',false);
medianMask=median(imMask,5);
datasetName=fullfile(CONSTANTS.imageData.imageDir,CONSTANTS.imageData.DatasetName);
end
function bwMembrane=getMembrane(CONSTANTS,channel,t)
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',channel,...
'timeRange',[t t], 'outType','single','prompt',false);
if isempty(im) || all(im(:)==0)
bwMembrane=[];
return;
end
imf=HIP.LoG(im,[3,3,0],[]);
imf(imf>0)=0;
imf=abs(imf);
bwMembrane=imbinarize(imf,adaptthresh(imf));
% bwMembrane=bwmorph(bwMembrane,'skel',Inf);
......@@ -3,7 +3,7 @@ function [functionHandle,params]=getSegmentHelper(conn,strAlgorithmType)
% set the default seg function and args
% default
if strcmp('segment',strAlgorithmType)
functionHandle=str2func('Segment.FrameSegment');
functionHandle=str2func('Segment.frameSegment');
else
functionHandle=[];
end
......
function bw=ljsThreshold(im,min_radius_pixels,min_area_pixels)
imOut=0*im;
for i=1:5
imf=HIP.Gaussian(im,[floor(size(im)/(128/2^(i)))+1,0],1,[]);
% imx=imf;
bw=(im - 1*(imf)>0);
imOut=imOut+bw;
end
se=strel('square',ceil(min_radius_pixels));
imo=imopen(imOut,se);
bw=logical(imo>0);
bw=bwareaopen(bw,min_area_pixels);
\ No newline at end of file
function bw=segReduce(bw,qLog,min_area_pixels,min_radius_pixels)
function bw=segReduce(bw,bwLog,segParams,min_area_pixels,min_radius_pixels)
[origL,num]=bwlabeln(bw);
% firstDiscovery is the gradient level at which a given connected component
% first finds a cell. We use this to constrain so that a very strong group
% of cells doesn't 'spawn' artifact cells from their intersection at lower
% gradient levels
firstDiscovery=zeros(1,num);
qMax=max(qLog(:));
% qMax=12;
bwErode=bw;
qL=[];
bwDilate=bwLog;
qL=origL;
% nDimension is 2 for 2-D, 3 for 3-D used for concat'ing
nDimension=length(size(bw));
for iq=qMax:-1:1
bwErode(qLog>=iq)=0;
ccErode=bwconncomp(bwErode);
cellSizes=cellfun(@length,ccErode.PixelIdxList);
% if is3D(bw)
% check that each CC achieves a distance of at least 0.5 * min_radius
% maybe bad idea -- holes in bwErode make distance unreliable?
% bwd=bwdist(~bwErode);
% rpdist=cellfun(@(x) max(bwd(x)),ccErode.PixelIdxList);
% idx=cellSizes<min_area_pixels | rpdist<1.5*min(min_radius_pixels);
% else
% rp=regionprops(ccErode,'MinorAxisLength');
% idx=cellSizes<min_area_pixels | [rp.MinorAxisLength]<2*min_radius_pixels;
% end
idx=cellSizes<min_area_pixels;
resetPixels=cell2mat(ccErode.PixelIdxList(idx)');
bwErode(resetPixels)=0;
if isempty(find(bwErode, 1))
if is3D(bw)
minKernelArea=ceil(min_area_pixels/2); %
else
minKernelArea=ceil(min_area_pixels/4); %
end
if is3D(bw)
se=strel('cube',3);
else
se=strel('disk',1);
end
for nDilate=1:2*min_radius_pixels
bwKernels=bw&~bwDilate;
bwKernels=bwareaopen(bwKernels,minKernelArea);
if isempty(find(bwKernels, 1))
break;
end
L=bwlabeln(bwErode);
L=bwlabeln(bwKernels);
qL=cat(nDimension+1,qL,L);
bwDilate=imdilate(bwDilate,se);
end
if nDimension==2
......@@ -60,25 +52,12 @@ for iq=niq:-1:1
rpdist=cellfun(@(x) max(bwd(x)),ccL.PixelIdxList);
% for each newL component, how many components from L are there?
for n=1:max(newL(:))
if iq>1 && (rpdist(n)<0.5*min(min_radius_pixels))
continue
end
% if iq>1 && (rpdist(n)<0.5*min(min_radius_pixels))
% continue
% end
idx=find(n==newL);
discoverL=unique(origL(idx));
discoverL(0==discoverL)=[];
if ~isempty(discoverL)
qDiscover=min(firstDiscovery(discoverL));
% if qDiscover-iq>2
% continue
% end
end
nc=numComponents(L,idx);
if nc<=1
if 0==nc
firstDiscovery(discoverL)=iq;
end
% keep this component from newL -- there are 0 or 1 children
bw(idx)=1;
else
......@@ -89,16 +68,9 @@ for iq=niq:-1:1
L=bwlabeln(bw);
end
4;
function nc=numComponents(L,idx)
labels=unique(L(idx));
labels(0==labels)=[];
nc=length(labels);
function b3D=is3D(im)
if length(size(im))==3
b3D=true;
else
b3D=false;
end
\ No newline at end of file
function [bw,bwIntensity,nLevels]=thresholdImages(im,imTexture,segParams,...
min_radius_pixels,min_area_pixels)
function [bw,bwLog]=thresholdImages(im,imLog,segParams, min_radius_pixels,min_area_pixels,...
medianMask)
bw=[];
nLevels=-1;
imLog(imLog<0)=0;
imLog=mat2gray(imLog);
bwIntensity=thresholdIntensity(im,segParams,min_radius_pixels,min_area_pixels);
% if no imTexture provided, just return bwIntensity
if isempty(imTexture)
if isfield(segParams,'sensitivity')
sensitivity=segParams.sensitivity;
else
sensitivity=0.5;
end
if is3D(im)
nsz=2*floor(size(im)/16)+1;
T=adaptthresh(imLog,0.5,'NeighborhoodSize',nsz,'statistic','gaussian');
bwLog=imbinarize(imLog,T);
bwLog=bwareaopen(bwLog,4*min_area_pixels);
nsz=4*floor(size(im)/16)+1;
T=adaptthresh(im,sensitivity,'NeighborhoodSize',nsz,'statistic','gaussian');
bw=imbinarize(im,T);
bw=bwareaopen(bw,min_area_pixels);
return;
end
alpha=segParams.alpha(1);
% else (not 3d)
nsz=4*floor(size(im)/16)+1;
T=adaptthresh(imLog,0.5,'NeighborhoodSize',nsz,'statistic','gaussian');
bwLog=imbinarize(imLog,T);
bwLog=bwareaopen(bwLog,4*min_area_pixels);
if segParams.isPhase
if segParams.brightLevels>=0
nLevels=segParams.brightLevels;
else
nLevels=2;
end
lm=alpha*multithresh(imTexture,nLevels);
q=imquantize(imTexture,lm);
bw=logical(q>nLevels);
else
if segParams.brightLevels>=0
n0=segParams.brightLevels;
else
n0=1;
end
if is3D(im)
ratioThresh=30;
se=strel('disk',1);
se2=strel('disk',ceil(min_radius_pixels));
if segParams.isPhase>0
T=adaptthresh(im,sensitivity,'foregroundpolarity','dark','statistic','gaussian');
else
ratioThresh=10;
T=adaptthresh(im,sensitivity,'foregroundpolarity','bright','statistic','gaussian');
end
for nLevels=n0:20
lm=alpha*multithresh(imTexture,nLevels);
q=imquantize(imTexture,lm);
bw=logical(q>nLevels);
if length(find(bw))/length(find(bwIntensity))<ratioThresh
break;
bw=imbinarize(im,T);
if segParams.isPhase>0
bw=imcomplement(bw);
end
bw=imclose(bw,se2);
if true==segParams.wellRadius
bwMask=imbinarize(medianMask,adaptthresh(medianMask,'statistic','gaussian'));
bw=bw&~bwMask;
bwh=imfill(bw,'holes');
nrep=0;
% if the mask forms a solid circle, the hole fill on bw will fill
% in the whole foreground. dilate mask until that doesn't happen...
while nrep<5 && length(find(bwh&~bw))/length(bw(:))>0.25
nrep=nrep+1;
bwMask=imdilate(bwMask,se);
bw=bw&~bwMask;
end
bw=imopen(bw,se);
bwLog=bwLog&~bwMask;
end
end
function bwIntensity=thresholdIntensity(im,segParams,min_radius_pixels,min_area_pixels)
if length(segParams.alpha)>1
alpha=segParams.alpha(2);
bw=imfill(bw,'holes');
else
alpha=segParams.alpha;
end
% first, threshold bwIntensity
if segParams.brightLevels>=0
nLevels=segParams.brightLevels;
else
nLevels=1;
% use 2x the default nhood size
T=adaptthresh(im,sensitivity,'NeighborhoodSize',4*floor(size(im)/16)+1,'statistic','gaussian');
bw=imbinarize(im,T);
end
bw=bwareaopen(bw,min_area_pixels);
if segParams.isPhase
bwIntensity=thresholdPhaseIntensity(im,alpha,min_radius_pixels,min_area_pixels);
return;
end
% else not phase
if alpha>1 && alpha==round(alpha)
lm=multithresh(im,alpha);
bwIntensity=imquantize(im,lm);
bwIntensity=logical(bwIntensity>1);
else
lm=alpha*multithresh(im,nLevels);
bwIntensity=imquantize(im,lm);
bwIntensity=logical(bwIntensity>nLevels);
end
function bwIntensity=thresholdPhaseIntensity(im,alpha,min_radius_pixels,min_area_pixels)
% dark centers
lm=alpha*multithresh(im,16);
bwIntensity=imquantize(im,lm);
bwIntensity=logical(bwIntensity<=2);
% bright blobs
lm=alpha*multithresh(im,3);
bwBright=imquantize(im,lm);
se=strel('disk',ceil(min_radius_pixels/2));
bwBright=logical(bwBright>=3);
bwBright=imerode(bwBright,se); % erode? or open?
bwIntensity=bwIntensity|bwBright;
%clean up and done
bwIntensity=bwIntensity|imfill(bwIntensity,'holes');
% bwIntensity=bwareaopen(bwIntensity,round(min_area_pixels/4));
function b3D=is3D(im)
if length(size(im))==3
b3D=true;
else
b3D=false;
end
......@@ -5,7 +5,7 @@ jsAlgorithmInfo=[];
jsAlgorithmInfo.name='2D+3D generic: texture+gradient partition';
jsAlgorithmInfo.type='segment';
jsAlgorithmInfo.params=segParams;
jsAlgorithmInfo.function='+Segment.FrameSegment_texture';
jsAlgorithmInfo.function='+Segment.frameSegment';
jsAlgorithmInfo.commandHost='matlab';
jsAlgorithmInfo=jsonencode(jsAlgorithmInfo);
tag=datestr(datetime());
......
......@@ -3,4 +3,4 @@ function ljsPath=getLeverjsPath(file)
thisFile=mfilename('fullpath');
[thisFolder,~,~]=fileparts(thisFile);
ljsPath=fullfile(thisFolder,'..\leverjs\',file);
\ No newline at end of file
ljsPath=fullfile(thisFolder,'..','leverjs',file);
\ No newline at end of file
......@@ -4,19 +4,17 @@
"type" : "segment",
"description" : "texture w/ recurvise partition",
"params" : {
"note":" - channels: use e.g. [3,2,1]. the first channel should be the primary -- the remainder contribute edges",
"channels":"[1]",
"channels":"1",
"membraneChannel":"-1",
"note_1":" - minimumRadius_um: for ensemble, use [start,step,end] e.g. [2,0.1,4]",
"minimumRadius_um":"4",
"useCuda":"false",
"note_alpha":" - smaller alpha values lower segmentation threshold",
"alpha":"1.0",
"note_sensi":"sensitivity is on [0,1] 0 least fg, 1 most fg",
"sensitivity":0.5,
"useCuda":"true",
"denoise": "true",
"note_NLM":"non-local means params (denoise==true): [h,search window,neighborhood,high-pass filter]",
"NLM":"[0.1,6,3,0]",
"NLM":"[0.1,12,1,0]",
"isPhase":false,
"note_brightLevels":" - brightLevels is the base number of foreground groups, e.g. {2} for {bg,noise,fg}. leave blank for default",
"brightLevels":"",
"wellRadius":-1,
"note_2":" - bCytoplasmic to true for objects with fine processes to preserve",
"bCytoplasmic":false,
......@@ -24,12 +22,9 @@
"nCores":4,
"storeEnsemble":"true",
"note_4":"bCircular adds a bit more weight to circular segmentations (ensemble only)",
"bCircular":"false",
"note_5":"clonalMotionComp sets tracking costs by clusters of objects moving together",
"clonalMotionComp":"false"
"bCircular":"false"
},
"function":"+Segment.FrameSegment_texture.m",
"function":"+Segment.frameSegment.m",
"commandHost":"matlab"
},
{
......@@ -96,6 +91,17 @@
},
"function":"+CellFeatures.netNuclearEnvelope.m",
"commandHost":"matlab"
},
{
"name": "membraneRatio",
"type" : "feature",
"description" : "% of fg membrane pixels along surface",
"params" : {
"name":"membraneRatio",
"channel":1,
"color":"#ff0000"
},
"function":"+CellFeatures.membraneRatio.m",
"commandHost":"matlab"
}
]
\ No newline at end of file
......@@ -6,17 +6,23 @@ if nargin<1
fprintf(1,'opening db from lever.state: %s\n',strDB);
fclose(fid);
end