Commit 1dfbe06b authored by Andrew Cohen's avatar Andrew Cohen

new matlab/CellFeatures for erk -- changed the calculation on inner/outer

parent 71b487eb
......@@ -6,6 +6,16 @@ params.color="#00ffff";
% check for param override
params=Segment.getParams(params,args);
% segmentation params -- for min radius
algs=Read.getAlgorithms(conn,'segment');
args=algs.params;
segParams=Segment.getDefaultSegParams();
segParams=Segment.getParams(segParams,args);
resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
min_radius_pixels = segParams.minimumRadius_um(1) ./ resolution_um;
min_radius_pixels=min_radius_pixels(1);
% todo -- handle anistropy in 3d!
cmd=['select cellID from tblCells where time=' num2str(t)];
q=ljsFetch(conn,cmd);
cellIDs=cell2mat(q);
......@@ -24,22 +34,28 @@ for i=1:length(cellIDs)
end
labelCells(idx)=c.cellID;
end
% find region surrounding segmentation
[bwd,idxd]=bwdist(labelCells);
idxd=labelCells(idxd);
idxd(bwd>5)=0;
dInterior=bwdist(~labelCells);
labelInner=labelCells;
[dExterior,idxd]=bwdist(labelCells);
labelOuter=labelCells(idxd);
% should this be 5 or min_radius_pixels? outside region?
labelOuter(dExterior>5)=0;
% should this be 0, or min_radius_pixels? i.e. buffer between?
labelOuter(dInterior>=min_radius_pixels)=0;
labelInner(dInterior<min_radius_pixels)=0;
for i=1:length(cellIDs)
if is3d
c=Read.getCell3D(conn,cellIDs(i));
idxSegmentation=sub2ind(size(im),c.pts(:,2),c.pts(:,1),c.pts(:,3));
else
c=Read.getCell(conn,cellIDs(i));
idxSegmentation=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
end
idxCyto=find(idxd==cellIDs(i));
ratio=getErkRatio(im,idxCyto,idxSegmentation);
idxInner=find(labelInner==c.cellID);
idxOuter=find(labelOuter==c.cellID);
ratio=getErkRatio(im,idxOuter,idxInner);
if isnan(ratio)
cmd=['delete from tblCellFeatures where cellID=' num2str(cellIDs(i))];
exec(conn,cmd);
......@@ -56,16 +72,16 @@ for i=1:length(cellIDs)
end
function ratio=getErkRatio(im,idxCyto,idxSegmentation)
function ratio=getErkRatio(im,idxOuter,idxInner)
ratio=NaN;
imCyto=im(idxCyto);
imCyto(0==imCyto)=[];
if length(imCyto)<5
imOuter=im(idxOuter);
imOuter(0==imOuter)=[];
if length(imOuter)<5
return
end
imSegmentation=im(idxSegmentation);
imSegmentation(0==imSegmentation)=[];
if length(imSegmentation)/length(idxSegmentation)<0.1
imInner=im(idxInner);
imInner(0==imInner)=[];
if length(imInner)/length(idxInner)<0.1
return
end
ratio=median(imCyto)/median(imSegmentation);
ratio=median(imInner)/median(imOuter);
......@@ -3,7 +3,8 @@ function Cell=getCell(conn,cellID)
Cell=[];
try
cmd=['SELECT cellID,time,trackID,centroid,u16pixels,maxRadius,channel,area,surface FROM tblCells WHERE cellID=' num2str(cellID)];
cmd=['SELECT cellID,time,trackID,centroid,u16pixels,maxRadius,channel,area,surface,segCC '...
' FROM tblCells WHERE cellID=' num2str(cellID)];
Q=ljsFetch(conn,cmd);
catch
Cell=Read.getCell3D(conn,cellID);
......@@ -25,5 +26,5 @@ Cell.maxRadius=Q{6};
Cell.channel=Q{7};
Cell.area=Q{8};
Cell.surface=jsondecode(Q{9});
Cell.segCC=Q{10};
......@@ -11,6 +11,7 @@ if nargin<3
CONSTANTS=Read.getConstants(conn);
end
segParams=Segment.getDefaultSegParams();
segParams=[];
segParams.draw=false;
segParams.minimumRadius_um=5;
......@@ -75,7 +76,7 @@ else
% 2D
if USE_CUDA
entropyRadius=round(min_radius_pixels);
kernel=HIP.MakeEllipsoidMask([entropyRadius,entropyRadius,0]);
kernel=HIP.MakeEllipsoidMask([entropyRadius,entropyRadius,0]);
imTexture=HIP.EntropyFilter(im,kernel,[]);
else
se=strel('square',2*floor(min_radius_pixels/2)+1);
......@@ -98,7 +99,7 @@ if DRAW
drawnow
end
[bw,bwIntensity,nLevels]=thresholdImages(im,imTexture,segParams);
[bw,bwIntensity,nLevels]=Segment.thresholdImages(im,imTexture,segParams);
fprintf(1,'frameSegmentTexture: min_radius_pixels=%s nLevels detected at %d\n',mat2str(min_radius_pixels,3),nLevels);
bw=bw|imfill(bw,'holes');
......@@ -141,7 +142,6 @@ if isempty(bw)
end
preShakeL=bwlabeln(bw);
bwdCells=bwdist(~bw);
[L,num]=Segment.allocateShake(bw,bw2,qLog,im);
[L2]=bwlabeln(bw2); % used to mark original cc for segCC field
preShakeL=L;
......@@ -202,14 +202,13 @@ for i=1:length(Cells)
end
Cells(i).nfg=length(find(bwIntensity(idx)));
Cells(i).ncc=length(find([Cells.segCC]==Cells(i).segCC));
Cells(i).maxd=max(bwdCells(idx));
if DRAW && ( (0==Cells(i).nfg && Cells(i).ncc>1)|| Cells(i).maxd<min_radius_pixels/2)
if DRAW && ( (0==Cells(i).nfg && Cells(i).ncc>1) )
plot(Cells(i).surface(:,1),Cells(i).surface(:,2),'color','k','linewidth',3);
end
end
if ~isempty(Cells)
idxParasite=find( (0==[Cells.nfg] & 1<[Cells.ncc]) | ([Cells.maxd]<min(min_radius_pixels)/2));
idxParasite=find( 0==[Cells.nfg] & 1<[Cells.ncc] );
Cells(idxParasite)=[];
end
......@@ -225,56 +224,6 @@ if DRAW
end
4;
function [bw,bwIntensity,nLevels]=thresholdImages(im,imTexture,segParams)
if length(segParams.alpha)>1
alpha=segParams.alpha(1);
alphaIntensity=segParams.alpha(2);
else
alpha=segParams.alpha;
alphaIntensity=segParams.alpha;
end
% first, threshold bwIntensity
if segParams.brightLevels>=0
nLevels=segParams.brightLevels;
else
nLevels=1;
end
lm=alphaIntensity*multithresh(im,nLevels);
bwIntensity=imquantize(im,lm);
bwIntensity=logical(bwIntensity>nLevels);
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;
else
ratioThresh=10;
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;
end
end
end
function Cells=getGradientDef(Cells,imLog,bPhase)
if bPhase
......
function segParams=getDefaultSegParams()
segParams=[];
segParams.draw=false;
segParams.minimumRadius_um=5;
segParams.channels=1; %[1CONSTANTS.imageData.NumberOfChannels];
segParams.isPhase=false;
segParams.wellRadius=-1;
segParams.useCuda=false;
segParams.nCores=4; % for ensemble seg only
segParams.storeEnsemble=true;
segParams.alpha=1.0;
segParams.denoise=true;
segParams.bCytoplasmic=false;
segParams.bCircular=false;
segParams.brightLevels=-1;
function [bw,bwIntensity,nLevels]=thresholdImages(im,imTexture,segParams)
bw=[];
bwIntensity=[];
nLevels=-1;
if length(segParams.alpha)>1
alpha=segParams.alpha(1);
alphaIntensity=segParams.alpha(2);
else
alpha=segParams.alpha;
alphaIntensity=segParams.alpha;
end
% first, threshold bwIntensity
if segParams.brightLevels>=0
nLevels=segParams.brightLevels;
else
nLevels=1;
end
lm=alphaIntensity*multithresh(im,nLevels);
bwIntensity=imquantize(im,lm);
bwIntensity=logical(bwIntensity>nLevels);
% if no imTexture provided, just return bwIntensity
if isempty(imTexture)
return;
end
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;
else
ratioThresh=10;
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;
end
end
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment