...
 
Commits (64)
...@@ -386,7 +386,7 @@ function getMCRpath() ...@@ -386,7 +386,7 @@ function getMCRpath()
return mcrPath; return mcrPath;
// unix 2017b installed // unix 2017b installed
mcrPath='/usr/local/MATLAB/R2017b/'; mcrPath='/usr/local/MATLAB/R2019b/';
if (fs.existsSync(mcrPath)) if (fs.existsSync(mcrPath))
return mcrPath; return mcrPath;
// unix 2017b MCR // unix 2017b MCR
......
...@@ -22,52 +22,6 @@ function anyFileInQ(processQ,leverFileList) ...@@ -22,52 +22,6 @@ function anyFileInQ(processQ,leverFileList)
return null; return null;
} // anyFileInQ } // 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) function batchSeg(leverFileList,fnCallback)
{ {
const process=require('./process.js'); const process=require('./process.js');
...@@ -80,16 +34,12 @@ function batchSeg(leverFileList,fnCallback) ...@@ -80,16 +34,12 @@ function batchSeg(leverFileList,fnCallback)
inUseFile=path.basename(inUseFile); inUseFile=path.basename(inUseFile);
fnCallback(409,"leverFile :: "+inUseFile+" :: is already in the process Q -- cannot batch segment"); fnCallback(409,"leverFile :: "+inUseFile+" :: is already in the process Q -- cannot batch segment");
return; return;
} }
let nComplete=0;
for (let i=0;i<leverFileList.length;i++) { for (let i=0;i<leverFileList.length;i++) {
var leverDB=new LeverDB(leverFileList[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:[]}; var leverCommand={leverFile:leverFileList[i],time:1,command:'batchSeg',params:[]};
process.addToQ(leverCommand); process.addToQ(leverCommand);
nComplete++;
if (nComplete>=leverFileList.length)
startBatchMonitor(leverFileList);
}); });
} }
fnCallback(200,"OK"); fnCallback(200,"OK");
......
...@@ -11,7 +11,7 @@ function optimize(dbFile,bQuiet=false) ...@@ -11,7 +11,7 @@ function optimize(dbFile,bQuiet=false)
var cmd; var cmd;
cmd='PRAGMA optimize'; cmd='PRAGMA optimize';
lDB.exec(cmd,function(){ lDB.exec(cmd,function(){
cmd='PRAGMA wal_checkpoint(full)'; cmd='PRAGMA journal_mode=WAL;PRAGMA wal_checkpoint(full)';
lDB.exec(cmd,function(){ lDB.exec(cmd,function(){
lDB.close(function (error) { lDB.close(function (error) {
if (null!=error) if (null!=error)
......
...@@ -3,7 +3,7 @@ const fs = require('fs'); ...@@ -3,7 +3,7 @@ const fs = require('fs');
// use sha1 so it will match our git commit tags // use sha1 so it will match our git commit tags
const hash = crypto.createHash('sha1'); 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) { if (3!=process.argv.length) {
console.error('USAGE: node goHash.js filename'); console.error('USAGE: node goHash.js filename');
......
...@@ -1297,8 +1297,16 @@ class LeverDB { ...@@ -1297,8 +1297,16 @@ class LeverDB {
fnCallback(vertices); 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 // read all the leverfiles in a folder, send one massive response
// unused, since it takes // unused, since it takes
// note - currently unused... // note - currently unused...
......
...@@ -701,6 +701,8 @@ function getLineageMap(rowsFamilies,rowsTracks,rowsPhenotypes) ...@@ -701,6 +701,8 @@ function getLineageMap(rowsFamilies,rowsTracks,rowsPhenotypes)
// the phenotypes are indexed into constants.phenotypes (TBD). phenotype 0 is death. // the phenotypes are indexed into constants.phenotypes (TBD). phenotype 0 is death.
for (i=0;i<rowsPhenotypes.length;i++) { for (i=0;i<rowsPhenotypes.length;i++) {
node=gLineageMap.get(rowsPhenotypes[i].trackID); node=gLineageMap.get(rowsPhenotypes[i].trackID);
if (undefined==node)
continue;
if (undefined==node.phenotypes) if (undefined==node.phenotypes)
node.phenotypes=[]; node.phenotypes=[];
if (undefined==node.phenotypes[rowsPhenotypes[i].time]) if (undefined==node.phenotypes[rowsPhenotypes[i].time])
......
...@@ -64,7 +64,7 @@ function checkToProcessQ() ...@@ -64,7 +64,7 @@ function checkToProcessQ()
internalAddToQ(leverCommand[0]); internalAddToQ(leverCommand[0]);
} }
} }
toProcessQtimer=setTimeout(checkToProcessQ,1000); toProcessQtimer=setTimeout(checkToProcessQ,2000);
} // checkToProcessQ } // checkToProcessQ
var toProcessQtimer; var toProcessQtimer;
...@@ -112,16 +112,11 @@ function internalAddToQ(leverCommand) ...@@ -112,16 +112,11 @@ function internalAddToQ(leverCommand)
// note - when batchSeg is active, we disallow subsequent batch segs // note - when batchSeg is active, we disallow subsequent batch segs
if (null!=batchSegFile && 'batchSeg'==leverCommand.command) { if (null!=batchSegFile && 'batchSeg'==leverCommand.command) {
addToProcessQ(leverCommand); addToProcessQ(leverCommand);
getLeverFileStatus(batchSegFile,function(leverFile,uiStatus){
if ('batchTrackPending'==uiStatus || 'idle'==uiStatus) {
batchSegFile=null;
}
});
return; return;
} }
if ('batchSeg'==leverCommand.command) { if ('batchSeg'==leverCommand.command) {
batchSegFile=leverCommand.leverFile; batchSegFile=leverCommand.leverFile;
} }
// first, see if the leverFile is in the processQ // first, see if the leverFile is in the processQ
for (var i=0;i<inProcessQ.length;i++) { for (var i=0;i<inProcessQ.length;i++) {
if (undefined===inProcessQ[i].leverFile || null===inProcessQ[i].leverFile) if (undefined===inProcessQ[i].leverFile || null===inProcessQ[i].leverFile)
...@@ -232,6 +227,7 @@ function getUserID(leverCommand,fnCallback) ...@@ -232,6 +227,7 @@ function getUserID(leverCommand,fnCallback)
if (undefined!==fnCallback) if (undefined!==fnCallback)
fnCallback(userID); fnCallback(userID);
} }
leverDB.close();
}); });
}); });
} }
...@@ -241,6 +237,7 @@ function getUserID(leverCommand,fnCallback) ...@@ -241,6 +237,7 @@ function getUserID(leverCommand,fnCallback)
userID=rows[0].userID; userID=rows[0].userID;
if (undefined!==fnCallback) if (undefined!==fnCallback)
fnCallback(userID); fnCallback(userID);
leverDB.close();
} }
}); });
...@@ -256,7 +253,6 @@ function addToQ(leverCommand) ...@@ -256,7 +253,6 @@ function addToQ(leverCommand)
if (gEditID>Number.MAX_SAFE_INTEGER) if (gEditID>Number.MAX_SAFE_INTEGER)
gEditID=1; gEditID=1;
} }
getUserID(leverCommand,function(userID){ getUserID(leverCommand,function(userID){
if (undefined===leverCommand.params) { if (undefined===leverCommand.params) {
leverCommand.params={}; leverCommand.params={};
...@@ -322,8 +318,13 @@ function addFiletoBatchSeg(leverCommand,nSlot) ...@@ -322,8 +318,13 @@ function addFiletoBatchSeg(leverCommand,nSlot)
}); });
inProcessQ[nSlot].childProc.stdout.on('data', (data) => { 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) => { inProcessQ[nSlot].childProc.stderr.on('data', (data) => {
...@@ -331,11 +332,21 @@ function addFiletoBatchSeg(leverCommand,nSlot) ...@@ -331,11 +332,21 @@ function addFiletoBatchSeg(leverCommand,nSlot)
}); });
inProcessQ[nSlot].childProc.on('close', (code) => { inProcessQ[nSlot].childProc.on('close', (code) => {
ljsLog.log('inProcessQ slot '+nSlot+' : exited with code ' + code); ljsLog.log('inProcessQ slot '+nSlot+' : exited with code ' + code);
var lf=inProcessQ[nSlot].leverFile;
inProcessQ[nSlot].childProc=null; inProcessQ[nSlot].childProc=null;
inProcessQ[nSlot].leverFile=null; inProcessQ[nSlot].leverFile=null;
inProcessQ[nSlot].status=""; 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 } // addFiletoBatchSeg
...@@ -344,7 +355,13 @@ function addFiletoBatchSeg(leverCommand,nSlot) ...@@ -344,7 +355,13 @@ function addFiletoBatchSeg(leverCommand,nSlot)
function addFileToProcess(leverCommand,nSlot) function addFileToProcess(leverCommand,nSlot)
{ {
if ("batchSeg"==leverCommand.command){ 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; return;
} }
const path=require('path'); const path=require('path');
...@@ -500,6 +517,7 @@ function removeID(editID) ...@@ -500,6 +517,7 @@ function removeID(editID)
function pauseOrKillAll(bKill,leverFile) function pauseOrKillAll(bKill,leverFile)
{ {
const spawn = require('child_process').spawn;
toProcessQ=[]; toProcessQ=[];
for (var i=0;i<inProcessQ.length;i++) { for (var i=0;i<inProcessQ.length;i++) {
...@@ -511,11 +529,21 @@ function pauseOrKillAll(bKill,leverFile) ...@@ -511,11 +529,21 @@ function pauseOrKillAll(bKill,leverFile)
if (bKill) { if (bKill) {
gCommandMap.delete(inProcessQ[i].ID); gCommandMap.delete(inProcessQ[i].ID);
if (null!==inProcessQ[i].childProc) if (null!==inProcessQ[i].childProc) {
inProcessQ[i].childProc.kill(); 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].childProc=null;
inProcessQ[i].leverFile=null; inProcessQ[i].leverFile=null;
batchSegFile=null; batchSegFile=null;
inProcessQ[nSlot].matlabPID=0;
} }
else { else {
// pause gets put right into Q // pause gets put right into Q
......
...@@ -950,12 +950,12 @@ function drawHulls(tHulls) ...@@ -950,12 +950,12 @@ function drawHulls(tHulls)
} }
// color by feature? // color by feature?
if (nFeature>=0) { if (nFeature>=0) {
var fmax=gFeatureExtrema[nFeature][1]; // var fmax=gFeatureExtrema[nFeature][1];
var fmin=gFeatureExtrema[nFeature][0]; // var fmin=gFeatureExtrema[nFeature][0];
// scale so colors go from [0.2,1] // // scale so colors go from [0.2,1]
// puts it on [0,1] // // puts it on [0,1]
var fval=1 - (tHulls[i][CONSTANTS.cellFeatures[0].name] - fmin)/((fmax-fmin)); var fval=tHulls[i][CONSTANTS.cellFeatures[0].name];
// now put it on [0.1,1] // // now put it on [0.1,1]
fval=0.9*fval+0.1; fval=0.9*fval+0.1;
color = getHTMLcolor([fval,fval,fval]); color = getHTMLcolor([fval,fval,fval]);
} }
......
...@@ -48,7 +48,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1 ...@@ -48,7 +48,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1
sp.NLM(4)=commandList{labindex+idxCommandList,'hpf'}; sp.NLM(4)=commandList{labindex+idxCommandList,'hpf'};
sp.alpha=commandList{labindex+idxCommandList,'alpha'}; sp.alpha=commandList{labindex+idxCommandList,'alpha'};
sp.ensemble_minimumRadius_um=true; sp.ensemble_minimumRadius_um=true;
nc=Segment.FrameSegment_texture('',1,CONSTANTS,sp); nc=Segment.frameSegment('',1,CONSTANTS,sp);
else else
nc=[]; nc=[];
end end
......
...@@ -27,6 +27,15 @@ else ...@@ -27,6 +27,15 @@ else
end end
segParams.draw=false; segParams.draw=false;
commandList=[]; 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 % make the command list
for t=1:CONSTANTS.imageData.NumberOfFrames for t=1:CONSTANTS.imageData.NumberOfFrames
nc=[repmat(t,length(rgRadius),1),rgRadius']; nc=[repmat(t,length(rgRadius),1),rgRadius'];
...@@ -68,7 +77,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1 ...@@ -68,7 +77,7 @@ for idp=0:ceil(size(commandList,1)/p.NumWorkers)-1
sp=segParams; sp=segParams;
sp.minimumRadius_um=commandList(labindex+idxCommandList,2); sp.minimumRadius_um=commandList(labindex+idxCommandList,2);
if isempty(intersect(t,tSeg)) && t>0 && t<=CONSTANTS.imageData.NumberOfFrames 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 else
if ~isempty(intersect(t,tSeg)) if ~isempty(intersect(t,tSeg))
fprintf(1,' found existing seg, t=%d, idxCommandList=%d\n',t,idxCommandList+labindex); fprintf(1,' found existing seg, t=%d, idxCommandList=%d\n',t,idxCommandList+labindex);
......
...@@ -19,11 +19,16 @@ for ff=1:length(flist) ...@@ -19,11 +19,16 @@ for ff=1:length(flist)
Batch.batchSegment(strDB,[],false); Batch.batchSegment(strDB,[],false);
end end
if ispc
nodecmd='node.exe ';
else
nodecmd='node ';
end
ljsPath=getLeverjsPath('resegAndTrackDataset.js'); ljsPath=getLeverjsPath('resegAndTrackDataset.js');
if bTrack if bTrack
parfor ff=1:length(flist) parfor ff=1:length(flist)
strDB=fullfile(flist(ff).folder,flist(ff).name); strDB=fullfile(flist(ff).folder,flist(ff).name);
system(['node.exe ' ljsPath ' "' ... system([nodecmd ' ' ljsPath ' "' ...
strDB '" resegAll 1 "{}" ' ]); strDB '" resegAll 1 "{}" ' ]);
end end
end end
......
...@@ -5,7 +5,7 @@ exec(conn,['ALTER TABLE tblCellFeatures ADD COLUMN ' featureName ' NUMBER']); ...@@ -5,7 +5,7 @@ exec(conn,['ALTER TABLE tblCellFeatures ADD COLUMN ' featureName ' NUMBER']);
nf.name=featureName; nf.name=featureName;
nf.color=color; nf.color=color;
if ~isfield(CONSTANTS,'cellFeatures') if ~isfield(CONSTANTS,'cellFeatures') || isempty(CONSTANTS.cellFeatures)
CONSTANTS.cellFeatures=nf; CONSTANTS.cellFeatures=nf;
else else
if ~any(cell2mat(strfind({CONSTANTS.cellFeatures.name},nf.name))) 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; ensembleDraw=false;
if segParams.draw==true if segParams.draw==true
...@@ -37,7 +37,7 @@ parfor i=1:length(rgRadius) ...@@ -37,7 +37,7 @@ parfor i=1:length(rgRadius)
% note as of matlab2019b, we pass [] for the conn (1st) argument as the % 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 % param is unused as long as CONSTANTS are provided, and even passing
% it around in a parfor generates a warning % 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; ensembleCells{i}=cx;
end end
segParams.draw=ensembleDraw; segParams.draw=ensembleDraw;
......
...@@ -9,12 +9,12 @@ nctime=tic(); ...@@ -9,12 +9,12 @@ nctime=tic();
% initialize the bucket with the cells from the maximal radius % initialize the bucket with the cells from the maximal radius
% segmentation % segmentation
idx=1; idx=1;
while isempty(ensembleCells{idx}) && idx<length(ensembleCells) while isempty(ensembleCells{idx}) && idx<=length(ensembleCells)
idx=idx+1; idx=idx+1;
end end
if idx==length(ensembleCells) if idx>length(ensembleCells)
% only zero or one non-empty bucket -- nothing to nest % only zero or one non-empty bucket -- nothing to nest
nestedCells=ensembleCells{idx}; % could be [], or just cells from smallest radius nestedCells={};
return; return;
end end
for i=1:length(ensembleCells{idx}) for i=1:length(ensembleCells{idx})
...@@ -30,6 +30,9 @@ end ...@@ -30,6 +30,9 @@ end
% the two buckets % the two buckets
for i=idx+1:length(ensembleCells) for i=idx+1:length(ensembleCells)
iCells=ensembleCells{i}; iCells=ensembleCells{i};
if isempty(iCells)
continue
end
p=gcp('nocreate'); p=gcp('nocreate');
if isempty(p) if isempty(p)
[iBucketDest,iCells]=getBucketDest(iCells,rgRadius,nestedCells,szIm,i); [iBucketDest,iCells]=getBucketDest(iCells,rgRadius,nestedCells,szIm,i);
......
...@@ -15,18 +15,26 @@ else ...@@ -15,18 +15,26 @@ else
se=strel('disk',3); se=strel('disk',3);
end end
imBG=imerode(im,se); 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) for i=1:length(Cells)
if is3D(im) if is3D(im)
verts=round(Cells(i).verts); verts=round(Cells(i).verts);
verts=[verts(:,2),verts(:,1),verts(:,3)];
verts=max(verts,1); verts=max(verts,1);
verts=min(verts,size(im)); 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 else
idx=sub2ind(size(im),Cells(i).surface(:,2),Cells(i).surface(:,1)); idx=sub2ind(size(im),Cells(i).surface(:,2),Cells(i).surface(:,1));
end end
imBGmu=mean(imBG(idx)); idx(isnan(idx))=[];
minFG=graythresh(im)*segParams.alpha(1); imBGmu=median(imBG(idx));
minFG=max(T(idx));
imBGeff=1-imBGmu/minFG; % imBGeff=1-imBGmu/minFG; %
imBGeff=min(imBGeff,1); 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 % area, convexArea, circularity on each cell
function newCell=setCellFeatures(newCell,RP,n,segParams) function newCell=setCellFeatures(newCell,RP,n,segParams)
......
...@@ -55,7 +55,7 @@ for t=1:CONSTANTS.imageData.NumberOfFrames ...@@ -55,7 +55,7 @@ for t=1:CONSTANTS.imageData.NumberOfFrames
if length(idx)<5 if length(idx)<5
continue continue
end end
newCell=Segment.FrameSegment_3D_create(idx,size(bwMeta),metaChannel,t); newCell=Segment.frameSegment_create_3D(idx,size(bwMeta),metaChannel,t);
metaCells=[metaCells newCell]; metaCells=[metaCells newCell];
end end
end end
......
...@@ -2,12 +2,11 @@ ...@@ -2,12 +2,11 @@
UISERVER='http://localhost:4444/' UISERVER='http://localhost:4444/'
movieFile='pos1_rotate_seg.mp4'; movieFile='9drotation.mp4';
v=VideoWriter(movieFile,'MPEG-4'); v=VideoWriter(movieFile,'MPEG-4');
v.FrameRate=10; v.FrameRate=10;
open(v) open(v)
ui=webread([UISERVER 'ui']); ui=webread([UISERVER 'ui']);
ui.sidebar='none'; ui.sidebar='none';
ui.webToolbar='none'; ui.webToolbar='none';
...@@ -21,19 +20,6 @@ webwrite([UISERVER 'ui'],options,ui); ...@@ -21,19 +20,6 @@ webwrite([UISERVER 'ui'],options,ui);
time=1; time=1;
webwrite([UISERVER 'time/' num2str(time)] ,''); 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']); view=webread([UISERVER 'view']);
options = weboptions('MediaType','application/json'); options = weboptions('MediaType','application/json');
theta=1/pi; theta=1/pi;
...@@ -42,11 +28,8 @@ 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]; % 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 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]; 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=ry;
% worldRot=worldRot';
view.worldRot=worldRot(:); view.worldRot=worldRot(:);
% view.zoom=0.5;
webwrite([UISERVER 'view'],options,view); webwrite([UISERVER 'view'],options,view);
bComplete=false; bComplete=false;
while ~bComplete while ~bComplete
...@@ -56,9 +39,9 @@ for theta=0:0.1:2*pi ...@@ -56,9 +39,9 @@ for theta=0:0.1:2*pi
end end
end end
im=webread([UISERVER 'screenCap']); im=webread([UISERVER 'screenCap']);
% imagesc(im) imagesc(im)
% axis image axis image
writeVideo(v,im); % writeVideo(v,im);
drawnow drawnow
fprintf(1,'theta=%f complete\n',theta); fprintf(1,'theta=%f complete\n',theta);
end end
......
...@@ -25,7 +25,7 @@ end ...@@ -25,7 +25,7 @@ end
srcTracks=cell2mat(q); srcTracks=cell2mat(q);
% tracks that start in this frame % tracks that start in this frame
cmd=['select cellID,time from tblCells where cellID=trackID and time>='... 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)']; ' and cellID not in (select cellID_dst from tblManualEdits)'];
q=ljsFetch(conn,cmd); q=ljsFetch(conn,cmd);
if isempty(q) if isempty(q)
......
...@@ -15,7 +15,7 @@ editQ=ljsFetch(conn,cmd); ...@@ -15,7 +15,7 @@ editQ=ljsFetch(conn,cmd);
if ~isempty(editQ) if ~isempty(editQ)
idxUndo = cellfun(@(x) strcmp(x,'undo'),editQ(:,5)); idxUndo = cellfun(@(x) strcmp(x,'undo'),editQ(:,5));
idxUndone = [editQ{idxUndo,1}]' idxUndone = [editQ{idxUndo,1}]';
for i=1:size(editQ,1) for i=1:size(editQ,1)
if idxUndo(i) if idxUndo(i)
continue; continue;
......
...@@ -84,7 +84,7 @@ end ...@@ -84,7 +84,7 @@ end
segParams=Read.getSegmentationParams(conn); segParams=Read.getSegmentationParams(conn);
ptsReplicate_xy=Segment.FrameSegment_texture_getSplitPixels(conn,CONSTANTS,... ptsReplicate_xy=Segment.frameSegment_getSplitPixels(conn,CONSTANTS,...
idxPixels, channel, time, segParams); idxPixels, channel, time, segParams);
try try
...@@ -106,7 +106,7 @@ for kpts=1:max(idx) ...@@ -106,7 +106,7 @@ for kpts=1:max(idx)
return; return;
end end
ptsIdx=find(idx==kpts); 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); channel,time,conn,segParams);
newCell.segCC=Q{6}; newCell.segCC=Q{6};
Cells=[Cells newCell]; Cells=[Cells newCell];
......
function newCell=FrameSegment_texture_create(idxPixels, imSize, channel, time, ... function newCell=frameSegment_create(idxPixels, imSize, channel, time, ...
conn, params) conn, params)
if 3==length(imSize) && imSize(3)>1 if 3==length(imSize) && imSize(3)>1
% 3D % 3D
newCell=Segment.FrameSegment_3D_create(idxPixels, imSize, channel, time); newCell=Segment.frameSegment_create_3D(idxPixels, imSize, channel, time);
return; return;
end end
......
function newCell=FrameSegment_3D_create(idxPixels, imSize, channel, time, ... function newCell=FrameSegment_create_3D(idxPixels, imSize, channel, time, ...
conn, params) conn, params)
bwn=zeros(imSize); bwn=zeros(imSize);
bwn(idxPixels)=1; bwn(idxPixels)=1;
PAD=5; PAD=5;
% padd with extra zeros around the outside so our mesh covers boundaries % padd with extra zeros around the outside so our mesh covers boundaries
bpad = zeros(size(bwn)+10); bpad = zeros(size(bwn)+10);
bpad(PAD+1:end-PAD,PAD+1:end-PAD,PAD+1:end-PAD)=bwn; bpad(PAD+1:end-PAD,PAD+1:end-PAD,PAD+1:end-PAD)=bwn;
bwn = bpad; bwn = bpad;
% scale=[1,1,1]; % scale=[1,1,1];
scale=max(imSize)./imSize/8; scale=max(imSize)./imSize/8;
scale=max(min(scale,0.5),.01); scale=max(min(scale,0.5),.01);
% %
resize=round(scale.*size(bwn)); resize=round(scale.*size(bwn));
bwResize=imresize3(bwn,resize); bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize)); [faces,verts]=isosurface(bwResize,graythresh(bwResize));
if size(verts,1)<10 if size(verts,1)<10
resize=round(2*scale.*size(bwn)); resize=round(2*scale.*size(bwn));
bwResize=imresize3(bwn,resize); bwResize=imresize3(bwn,resize);
[faces,verts]=isosurface(bwResize,graythresh(bwResize)); [faces,verts]=isosurface(bwResize,graythresh(bwResize));
end end
norms = isonormals(bwResize,verts); norms = isonormals(bwResize,verts);
% rescale verts back to % rescale verts back to
unscale=(size(bwn)./size(bwResize)); unscale=(size(bwn)./size(bwResize));
% put unscale on x,y instead of r,c % put unscale on x,y instead of r,c
unscale=[unscale(2),unscale(1),unscale(3)]; unscale=[unscale(2),unscale(1),unscale(3)];
verts=(verts).*unscale-0.5*unscale; verts=(verts).*unscale-0.5*unscale;
edges=Segment.MakeEdges(faces); edges=Segment.MakeEdges(faces);
% make edges and faces zero indexed % make edges and faces zero indexed
edges=edges-1; edges=edges-1;
faces=faces-1; faces=faces-1;
% subtract extract for padding % subtract extract for padding
verts=verts-PAD; verts=verts-PAD;
maxRad = verts-repmat(mean(verts,1),size(verts,1),1); maxRad = verts-repmat(mean(verts,1),size(verts,1),1);
maxRad=sum(abs(maxRad),2); maxRad=sum(abs(maxRad),2);
maxRad=max(maxRad); maxRad=max(maxRad);
newCell=[]; newCell=[];
newCell.time=time; newCell.time=time;
newCell.centroid=[mean(verts,1)]-1; newCell.centroid=[mean(verts,1)]-1;
newCell.edges=edges; newCell.edges=edges;
newCell.faces=faces; newCell.faces=faces;
newCell.verts=verts; newCell.verts=verts;
newCell.normals=norms; newCell.normals=norms;
xyzPts=[]; xyzPts=[];
% note that xyzPts are on the correct range, as idxPixels is unpadded... % note that xyzPts are on the correct range, as idxPixels is unpadded...
[xyzPts(:,2),xyzPts(:,1),xyzPts(:,3)]=ind2sub(imSize,idxPixels); [xyzPts(:,2),xyzPts(:,1),xyzPts(:,3)]=ind2sub(imSize,idxPixels);
newCell.pts=uint16(xyzPts); newCell.pts=uint16(xyzPts);
newCell.maxRadius=maxRad; newCell.maxRadius=maxRad;
newCell.channel=channel; 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) params)
segParams=Segment.getDefaultSegParams(); segParams=Segment.getDefaultSegParams();
......
...@@ -2,19 +2,16 @@ function segParams=getDefaultSegParams() ...@@ -2,19 +2,16 @@ function segParams=getDefaultSegParams()
segParams=[]; segParams=[];
segParams.draw=false; segParams.draw=false;
segParams.minimumRadius_um=[1.5,0.25,2.5]; 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.isPhase=false;
segParams.wellRadius=-1; segParams.wellRadius=-1;
segParams.useCuda=true; segParams.useCuda=true;
segParams.nCores=feature('numcores'); % for ensemble seg only segParams.nCores=feature('numcores'); % for ensemble seg only
segParams.storeEnsemble=true; segParams.storeEnsemble=true;
segParams.alpha=1.0; segParams.sensitivity=0.5;
segParams.denoise=true; segParams.denoise=true;
segParams.bCytoplasmic=false; segParams.bCytoplasmic=false;
segParams.bCircular=false; segParams.bCircular=false;
segParams.brightLevels=-1;
segParams.clonalMotionComp=false;
% advanced params
%
% NLM: h,search_window_radius,neighborhood_radius,high_pass_filter % NLM: h,search_window_radius,neighborhood_radius,high_pass_filter
segParams.NLM=[0.1,12,1,false]; 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') if ~exist('bFilter','var')
bFilter=false; bFilter=false;
...@@ -8,17 +8,22 @@ resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel ...@@ -8,17 +8,22 @@ resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
min_radius_pixels = segParams.minimumRadius_um ./ resolution_um; min_radius_pixels = segParams.minimumRadius_um ./ resolution_um;
global USE_CUDA global USE_CUDA
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',channel,'timeRange',[t t], 'outType','single','prompt',false); 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) if isempty(im)
qLog=[]; imLog=[];
imLogRaw=[];
return return
end end
imRaw=im;
im=denoise(im,segParams); im=denoise(im,segParams);
if 1==nargout if 1==nargout
return return
end end
if is3D(im) if is3D(im)
% 3D % 3D
if 1==nargout if 1==nargout
...@@ -26,43 +31,42 @@ if is3D(im) ...@@ -26,43 +31,42 @@ if is3D(im)
end end
if USE_CUDA if USE_CUDA
logRadius=min_radius_pixels.*(1/sqrt(3)); logRadius=min_radius_pixels.*(1/sqrt(3));
imLog=HIP.LoG(im,logRadius,[]); imLog=HIP.LoG(imRaw,logRadius,[]);
else else
% use difference of gaussian approximation % use difference of gaussian approximation
logRadius=min_radius_pixels; logRadius=min_radius_pixels;
imd1=imgaussfilt3(im,sqrt(2).*logRadius); imd1=imgaussfilt3(imRaw,sqrt(2).*logRadius);
imd2=imgaussfilt3(im,logRadius./sqrt(2)); imd2=imgaussfilt3(imRaw,logRadius./sqrt(2));
imLog=imd1-imd2; imLog=imd1-imd2;
end end
imLogRaw=imLog;
else else
% 2D % 2D
if USE_CUDA if USE_CUDA
logRadius=(1/sqrt(2)).*min_radius_pixels; logRadius=(1/sqrt(2)).*min_radius_pixels;
logRadius(2)=logRadius(1);
logRadius(3)=0; logRadius(3)=0;
imLog=HIP.LoG(im,logRadius,[]); imLog=HIP.LoG(imRaw,logRadius,[]);
else else
if length(min_radius_pixels)>1 if length(min_radius_pixels)>1
logRadius=0.5*min(min_radius_pixels(1:2)); logRadius=0.5*min(min_radius_pixels(1:2));
else else
logRadius=0.5*min_radius_pixels ; logRadius=0.5*min_radius_pixels ;
end end
szFilter=round(size(im)/3); szFilter=round(size(imRaw)/3);
h=fspecial('log',szFilter,logRadius); h=fspecial('log',szFilter,logRadius);
imLog=imfilter(im,h,'replicate'); imLog=imfilter(imRaw,h,'replicate');
end end
end end
imLogRaw=imLog;
imLog=mat2gray(imLog); % note - DO NOT mat2gray imLog here. Later processing can use the negative
lm=multithresh(imLog,15); % or positive portions for inside/outside info. Once we mat2gray that's
qLog=imquantize(imLog,lm); % gone.
if segParams.isPhase % LEFT HERE AS CAUTIONARY TALE imLog=mat2gray(imLog);
qLog=max(qLog(:))+1-qLog; %
end 4;
function im=denoise(im,segParams) function im=denoise(im,segParams)
if false==segParams.denoise if false==segParams.denoise
im=mat2gray(im); im=mat2gray(im);
return return
...@@ -88,7 +92,7 @@ if is3D(im) ...@@ -88,7 +92,7 @@ if is3D(im)
imx=imx-HIP.Cuda.Gaussian(im,szFilter,1,[]); imx=imx-HIP.Cuda.Gaussian(im,szFilter,1,[]);
end end
else else
imx = medfilt3(im)-imgaussfilt3(im,szFilter); imx = medfilt3(im);
end end
else else
% 2D % 2D
...@@ -98,10 +102,10 @@ else ...@@ -98,10 +102,10 @@ else
imx=imx-HIP.Cuda.Gaussian(im,[szFilter,1],1,[]); imx=imx-HIP.Cuda.Gaussian(im,[szFilter,1],1,[]);
end end
else else
imx=medfilt2(im)-imgaussfilt(im,szFilter(1:2)); imx=imnlmfilt(im);
end end
end end
imx=max(imx,0);
im=mat2gray(imx); 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) ...@@ -3,7 +3,7 @@ function [functionHandle,params]=getSegmentHelper(conn,strAlgorithmType)
% set the default seg function and args % set the default seg function and args
% default % default
if strcmp('segment',strAlgorithmType) if strcmp('segment',strAlgorithmType)
functionHandle=str2func('Segment.FrameSegment'); functionHandle=str2func('Segment.frameSegment');
else else
functionHandle=[]; functionHandle=[];
end 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); [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(:)); bwDilate=bwLog;
% qMax=12; qL=origL;
bwErode=bw;
qL=[];
% nDimension is 2 for 2-D, 3 for 3-D used for concat'ing % nDimension is 2 for 2-D, 3 for 3-D used for concat'ing
nDimension=length(size(bw)); nDimension=length(size(bw));
for iq=qMax:-1:1 if is3D(bw)
bwErode(qLog>=iq)=0; minKernelArea=ceil(min_area_pixels/2); %
ccErode=bwconncomp(bwErode); else
cellSizes=cellfun(@length,ccErode.PixelIdxList); minKernelArea=ceil(min_area_pixels/4); %
% if is3D(bw) end
% check that each CC achieves a distance of at least 0.5 * min_radius
% maybe bad idea -- holes in bwErode make distance unreliable? if is3D(bw)
% bwd=bwdist(~bwErode); se=strel('cube',3);
% rpdist=cellfun(@(x) max(bwd(x)),ccErode.PixelIdxList); else
% idx=cellSizes<min_area_pixels | rpdist<1.5*min(min_radius_pixels); se=strel('disk',1);
% else end
% rp=regionprops(ccErode,'MinorAxisLength');
% idx=cellSizes<min_area_pixels | [rp.MinorAxisLength]<2*min_radius_pixels; for nDilate=1:2*min_radius_pixels
% end bwKernels=bw&~bwDilate;
idx=cellSizes<min_area_pixels; bwKernels=bwareaopen(bwKernels,minKernelArea);
resetPixels=cell2mat(ccErode.PixelIdxList(idx)'); if isempty(find(bwKernels, 1))
bwErode(resetPixels)=0;
if isempty(find(bwErode, 1))
break; break;
end end
L=bwlabeln(bwErode); L=bwlabeln(bwKernels);
qL=cat(nDimension+1,qL,L); qL=cat(nDimension+1,qL,L);
bwDilate=imdilate(bwDilate,se);
end end
if nDimension==2 if nDimension==2
...@@ -60,25 +52,12 @@ for iq=niq:-1:1 ...@@ -60,25 +52,12 @@ for iq=niq:-1:1
rpdist=cellfun(@(x) max(bwd(x)),ccL.PixelIdxList); rpdist=cellfun(@(x) max(bwd(x)),ccL.PixelIdxList);
% for each newL component, how many components from L are there? % for each newL component, how many components from L are there?
for n=1:max(newL(:)) for n=1:max(newL(:))
if iq>1 && (rpdist(n)<0.5*min(min_radius_pixels)) % if iq>1 && (rpdist(n)<0.5*min(min_radius_pixels))
continue % continue
end % end
idx=find(n==newL); 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); nc=numComponents(L,idx);
if nc<=1 if nc<=1
if 0==nc
firstDiscovery(discoverL)=iq;
end
% keep this component from newL -- there are 0 or 1 children % keep this component from newL -- there are 0 or 1 children
bw(idx)=1; bw(idx)=1;
else else
...@@ -89,16 +68,9 @@ for iq=niq:-1:1 ...@@ -89,16 +68,9 @@ for iq=niq:-1:1
L=bwlabeln(bw); L=bwlabeln(bw);
end end
4;
function nc=numComponents(L,idx) function nc=numComponents(L,idx)
labels=unique(L(idx)); labels=unique(L(idx));
labels(0==labels)=[]; labels(0==labels)=[];
nc=length(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,... function [bw,bwLog]=thresholdImages(im,imLog,segParams, min_radius_pixels,min_area_pixels,...
min_radius_pixels,min_area_pixels) medianMask)
bw=[]; imLog(imLog<0)=0;
nLevels=-1; imLog=mat2gray(imLog);
bwIntensity=thresholdIntensity(im,segParams,min_radius_pixels,min_area_pixels); if isfield(segParams,'sensitivity')
% if no imTexture provided, just return bwIntensity sensitivity=segParams.sensitivity;
if isempty(imTexture) 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; return;
end 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.isPhase
if segParams.brightLevels>=0 se=strel('disk',1);
nLevels=segParams.brightLevels; se2=strel('disk',ceil(min_radius_pixels));
else if segParams.isPhase>0
nLevels=2; T=adaptthresh(im,sensitivity,'foregroundpolarity','dark','statistic','gaussian');
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;
else else
ratioThresh=10; T=adaptthresh(im,sensitivity,'foregroundpolarity','bright','statistic','gaussian');
end end
for nLevels=n0:20 bw=imbinarize(im,T);
lm=alpha*multithresh(imTexture,nLevels); if segParams.isPhase>0
q=imquantize(imTexture,lm); bw=imcomplement(bw);
bw=logical(q>nLevels); end
if length(find(bw))/length(find(bwIntensity))<ratioThresh bw=imclose(bw,se2);
break; 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 end
bw=imopen(bw,se);
bwLog=bwLog&~bwMask;
end end
end bw=imfill(bw,'holes');
function bwIntensity=thresholdIntensity(im,segParams,min_radius_pixels,min_area_pixels)
if length(segParams.alpha)>1
alpha=segParams.alpha(2);
else else
alpha=segParams.alpha; % use 2x the default nhood size
end T=adaptthresh(im,sensitivity,'NeighborhoodSize',4*floor(size(im)/16)+1,'statistic','gaussian');
% first, threshold bwIntensity bw=imbinarize(im,T);
if segParams.brightLevels>=0
nLevels=segParams.brightLevels;
else
nLevels=1;
end