Commit 671ca2b7 authored by actb's avatar actb

intensity based thresholding

parent 635ca11d
......@@ -25,6 +25,7 @@ for i=1:length(Cells)
else
idx=sub2ind(size(im),Cells(i).surface(:,2),Cells(i).surface(:,1));
end
idx(isnan(idx))=[];
imBGmu=mean(imBG(idx));
minFG=graythresh(im)*segParams.alpha(1);
imBGeff=1-imBGmu/minFG; %
......
......@@ -35,17 +35,12 @@ 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,medianMask);
[im,qLog,bwLog]=Segment.getImages(CONSTANTS,t,chan,segParams,medianMask);
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,medianMask);
qLog=qLog+channel_qLog;
end
if is3D(im)
min_area_pixels = 4/3*pi* mean(min_radius_pixels)^3;
......@@ -55,52 +50,20 @@ else
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);axis image;brighten(0.5);hold on
drawnow
end
[bw,bwIntensity,nLevels]=Segment.thresholdImages(im,imTexture,segParams,...
min_radius_pixels,min_area_pixels,medianMask);
sz=sprintf('frameSegmentTexture: min_radius_pixels=%s nLevels detected at %d\n',mat2str(min_radius_pixels,3),nLevels);
ljsLog(sz,3);
if ~segParams.isPhase
bw=bw|imfill(bw,'holes');
end
bw=imbinarize(im,adaptthresh(im,0.5,'statistic','gaussian'));
if segParams.wellRadius>1
% mask all points at r>radius
center=size(im)/2;
......@@ -111,62 +74,33 @@ if segParams.wellRadius>1
bwMask=logical(rMask<segParams.wellRadius);
bw(~bwMask)=0;
end
%
% if ~segParams.bCytoplasmic
% bw=removeOuterRadius(bw,CONSTANTS,min_radius_pixels);
% 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
% 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=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;
bw=bwareaopen(bw,min_area_pixels);
bw2=bw;
bw(bwLog)=0;
if ~segParams.isPhase
bw=bw|imfill(bw,'holes');
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);
[d,dNN]=bwdist(bw);
[L,num]=bwlabel(bw);
bwAssign=bw2&~bw;
bwAssign(d>min_radius_pixels)=0;
idxAssign=find(bwAssign);
L(idxAssign)=L(dNN(idxAssign));
num=max(L(:));
% tic
[L,num]=Segment.allocateShake(bw,bw2,qLog,min_radius_pixels);
% tShake=toc;
% fprintf(1,'allocateShake: t=%f\n',tShake);
% [L,num]=bwlabeln(bw);
% [L,num]=Segment.allocateShake(bw,bw2,qLog,min_radius_pixels);
[L2]=bwlabeln(bw2); % used to mark original cc for segCC field
preShakeL=L;
preShakeL(~bw)=0;
if bEnsemble
if size(bw,3)>1
......@@ -191,7 +125,6 @@ for n=1:num
segL=segL(1);
end
newCell.segCC=segL;
newCell.idxPreShake=find(preShakeL==n);
if bEnsemble
newCell=Ensemble.setCellFeatures(newCell,rp,n,segParams);
end
......@@ -208,41 +141,41 @@ for n=1:num
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(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
% 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
% % 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);
......@@ -255,7 +188,7 @@ if bEnsemble
end
end
end
end
tElapsed=toc;
sz=sprintf('segmented frame %d found %d hulls, useCuda=%d, time=%f\n',t,length(Cells),USE_CUDA,tElapsed);
......
function [im,qLog,imLogRaw]=getImages(CONSTANTS,t,channel,segParams,medianMask)
function [im,qLog,bwLog]=getImages(CONSTANTS,t,channel,segParams,medianMask)
if ~exist('bFilter','var')
bFilter=false;
......@@ -40,7 +40,6 @@ if is3D(im)
imd2=imgaussfilt3(im,logRadius./sqrt(2));
imLog=imd1-imd2;
end
imLogRaw=imLog;
else
% 2D
if USE_CUDA
......@@ -60,15 +59,18 @@ else
end
imLogRaw=imLog;
imLog=mat2gray(imLog);
%
lm=multithresh(imLog,15);
qLog=imquantize(imLog,lm);
if segParams.isPhase<0
% dark center phase -- need log to be small inside, large outside, so
% reverse
qLog=max(qLog(:))+1-qLog;
end
min_radius_pixels=min_radius_pixels(1);
min_area_pixels = min_radius_pixels^2 * pi;
min_area_pixels = max(1,round(min_area_pixels));
imLog(imLog<median(imLog(:)))=0;
imLog=mat2gray(imLog);
bwLog=imbinarize(imLog,adaptthresh(imLog,0.5,'statistic','gaussian'));
bwLog=bwareaopen(bwLog,2*min_area_pixels);
se=strel('disk',ceil(min_radius_pixels(1)));
bwLog=imclose(bwLog,se);
function im=denoise(im,segParams)
......
function [bw,bwIntensity,nLevels]=thresholdImages(im,imTexture,segParams,...
function [bw]=thresholdImages(im,imTexture,segParams,...
min_radius_pixels,min_area_pixels,medianMask)
bw=[];
......
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