thresholdImages : cleanLoG to denoise LoG

matlab/+Segment/thresholdImages.m

 
-function [bw,bwLog]=thresholdImages(im,imLog,segParams, min_radius_pixels,min_area_pixels,...
+function [bw,bwLoG]=thresholdImages(im,imLoG,segParams, min_radius_pixels,min_area_pixels,...
     medianMask)
 
-imLog(imLog<0)=0;
-imLog=mat2gray(imLog);
+imLoG(imLoG<0)=0;
+imLoG=mat2gray(imLoG);
 
 if isfield(segParams,'sensitivity')
     sensitivity=segParams.sensitivity;
@@ -12,9 +12,9 @@ else
 end
 
 if sensitivity>1 && round(sensitivity)==sensitivity
-    lm=multithresh(imLog,1);
-    bwLog=logical(imLog>lm(end));
-    bwLog=Segment.aiAreaOpen(bwLog,min_radius_pixels);
+    lm=multithresh(imLoG,1);
+    bwLoG=logical(imLoG>lm(end));
+    bwLoG=Segment.aiAreaOpen(bwLoG,min_radius_pixels);
     lm=multithresh(im,sensitivity);
     bw=logical(im>lm(end));
     bw=Segment.aiAreaOpen(bw,min_radius_pixels);
@@ -22,11 +22,11 @@ if sensitivity>1 && round(sensitivity)==sensitivity
 end
 
 nsz=4*floor(size(im)/16)+1;
-T=adaptthresh(imLog,0.5,'NeighborhoodSize',nsz,'statistic','gaussian');
-bwLog=imbinarize(imLog,T);
+T=adaptthresh(imLoG,0.5,'NeighborhoodSize',nsz,'statistic','gaussian');
+bwLoG=imbinarize(imLoG,T);
 
 if segParams.isPhase
-    [bw,bwLog]=phaseThreshold(im,bwLog,min_radius_pixels,segParams,...
+    [bw,bwLoG]=phaseThreshold(im,bwLoG,min_radius_pixels,segParams,...
         sensitivity,medianMask);
 else
     if is3D(im)
@@ -39,11 +39,13 @@ else
     bw=imbinarize(im,T);
    
 end
-bwLog=Segment.aiAreaOpen(bwLog,2*min_radius_pixels,true);
+% bwLog=Segment.aiAreaOpen(bwLog,2*min_radius_pixels,true);
+% bwLog=bwareaopen(bwLog,floor(2*pi*min_radius_pixels));
+bwLoG=cleanLoG(bwLoG);
 bw=Segment.aiAreaOpen(bw,min_radius_pixels);
 4;
 
-function [bw,bwLog]=phaseThreshold(im,bwLog,min_radius_pixels,segParams,...
+function [bw,bwLoG]=phaseThreshold(im,bwLoG,min_radius_pixels,segParams,...
     sensitivity,medianMask)
 se=strel('disk',1);
 rClose=max(1,ceil(min_radius_pixels/2));
@@ -59,19 +61,19 @@ if segParams.isPhase>0
 end
 
 if true==segParams.wellRadius
-    [bw,bwLog]=thresholdMask(medianMask,bw,bwLog,se,se2);
+    [bw,bwLoG]=thresholdMask(medianMask,bw,bwLoG,se,se2);
 end
 % fill in holes
 bw=imclose(bw,se2);
 bw=imopen(bw,se);
 bw=imfill(bw,'holes');
 
-bwLog=imclose(bwLog,se2);
+bwLoG=imclose(bwLoG,se2);
 4;
 
 
 
-function [bw,bwLog]=thresholdMask(medianMask,bw,bwLog,se,se2)
+function [bw,bwLoG]=thresholdMask(medianMask,bw,bwLoG,se,se2)
 
 imf=stdfilt(medianMask,true(5));
 bwMask=imbinarize(imf);
@@ -90,8 +92,14 @@ while nrep<5 && length(find(bwh&~bw))/length(bw(:))>0.25
 end
 bw=imfill(bw,'holes');
 bw=imopen(bw,se);
-bwLog=bwLog&~bwMask;
-
+bwLoG=bwLoG&~bwMask;
 
+% remove spurious LoG noise
+function bwLoG=cleanLoG(bwLoG)
 
+rp=regionprops(bwLoG,'ConvexArea','Area','PixelIdxList','filledarea');
+convDef=[rp.ConvexArea]./[rp.Area];
+idxReset=find(convDef<1.1);
+pixelIdxReset=vertcat(rp(idxReset).PixelIdxList);
+bwLoG(pixelIdxReset)=0;