superSeggerOpti generates the initial segmentation of rod-shaped cells.
It uses a local minimum filter (similar to a median filter) to enhance
contrast and then uses Matlab's WATERSHED command to generate
cell boundaries. The spurious boundaries (i.e., those that lie in the
cell interiors) are removed by an intensity thresholding routine
on each boundary. Any real boundaries incorrectly removed
by this thresholding are added back by an iterative algorithm that
uses knowledge of cell shape to determine which regions are missing
boundaries.
INPUT :
phaseOrData : phase image or data file to be used
mask : cell mask, given externally or calculated with band-pass filter
disp_flag : display flag
CONST : segmentation constants
adapt_flag : break up regions that are too big to be cells
header : string displayed with infromation
crop_box : information about alignement of the image
OUTPUT :
data.segs : defined below
data.mask_bg : a binary image in which all background (non-cell) pixels are masked
data.mask_cell : cell mask, a binary image the same size as phase in
which each cell is masked by a connected region of white pixels
data.phase : Original phase image
A : scoring vector optimized for different cells and imaging conditions
segs.
phaseMagic: % phase image processed with magicContrast only
segs_good: % on segments, image of the boundaries between cells that the program
has determined are correct (i.e., not spurious).
segs_bad: % off segments, image of program-determined spurious boundaries between cells
segs_3n: % an image of all of boundary intersections, segments that cannot be switched off
info: % segment parameters that are used to generate the raw
score, looke below
segs_label: % bwlabel of good and bad segs.
score: % cell scores for regions
scoreRaw: % raw scores for segments
props: % segement properties for segments
seg.info(:,1) : the minimum phase intensity on the seg
seg.info(:,2) : the mean phase intensity on the seg
seg.info(:,3) : area of the seg
seg.info(:,4) : the mean second d of the phase normal to the seg
seg.info(:,5) : second d of the phase normal to the seg at the min pixel
seg.info(:,6) : second d of the phase parallel to the seg at the min pixel
seg.info(:,7) and seg_info(:,8) : min and max area of neighboring regions
seg.info(:,9) and seg_info(:,10) : min and max lengths of the minor axis of the neighboring regions
seg.info(:,11) and seg_info(:,12) : min and max lengths of the major axis of the neighboring regions
seg.info(:,11) : length of minor axis
seg.info(:,12) : length of major axis
seg.info(:,13) : square of length of major axis
seg.info(:,16) : max length of region projected onto the major axis
segment
seg.info(:,17) : min length of region projected onto the major axis
segment
seg.info(:,18) : max length of region projected onto the minor axis
segment
seg.info(:,19) : min length of region projected onto the minor axis
segment
The output images are related by
mask_cell = mask_bg .* (~segs_good) .* (~segs_3n);
Copyright (C) 2016 Wiggins Lab
Written by Stella Stylianidou & Paul Wiggins.
University of Washington, 2016
This file is part of SuperSegger.
SuperSegger is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
SuperSegger is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with SuperSegger. If not, see <http://www.gnu.org/licenses/>.0001 function [data,A] = superSeggerOpti(phaseOrData, mask, disp_flag, CONST, adapt_flag, header, crop_box) 0002 % superSeggerOpti generates the initial segmentation of rod-shaped cells. 0003 % It uses a local minimum filter (similar to a median filter) to enhance 0004 % contrast and then uses Matlab's WATERSHED command to generate 0005 % cell boundaries. The spurious boundaries (i.e., those that lie in the 0006 % cell interiors) are removed by an intensity thresholding routine 0007 % on each boundary. Any real boundaries incorrectly removed 0008 % by this thresholding are added back by an iterative algorithm that 0009 % uses knowledge of cell shape to determine which regions are missing 0010 % boundaries. 0011 % 0012 % INPUT : 0013 % phaseOrData : phase image or data file to be used 0014 % mask : cell mask, given externally or calculated with band-pass filter 0015 % disp_flag : display flag 0016 % CONST : segmentation constants 0017 % adapt_flag : break up regions that are too big to be cells 0018 % header : string displayed with infromation 0019 % crop_box : information about alignement of the image 0020 % 0021 % OUTPUT : 0022 % data.segs : defined below 0023 % data.mask_bg : a binary image in which all background (non-cell) pixels are masked 0024 % data.mask_cell : cell mask, a binary image the same size as phase in 0025 % which each cell is masked by a connected region of white pixels 0026 % data.phase : Original phase image 0027 % A : scoring vector optimized for different cells and imaging conditions 0028 % 0029 % segs. 0030 % phaseMagic: % phase image processed with magicContrast only 0031 % segs_good: % on segments, image of the boundaries between cells that the program 0032 % has determined are correct (i.e., not spurious). 0033 % segs_bad: % off segments, image of program-determined spurious boundaries between cells 0034 % segs_3n: % an image of all of boundary intersections, segments that cannot be switched off 0035 % info: % segment parameters that are used to generate the raw 0036 % score, looke below 0037 % segs_label: % bwlabel of good and bad segs. 0038 % score: % cell scores for regions 0039 % scoreRaw: % raw scores for segments 0040 % props: % segement properties for segments 0041 % 0042 % 0043 % seg.info(:,1) : the minimum phase intensity on the seg 0044 % seg.info(:,2) : the mean phase intensity on the seg 0045 % seg.info(:,3) : area of the seg 0046 % seg.info(:,4) : the mean second d of the phase normal to the seg 0047 % seg.info(:,5) : second d of the phase normal to the seg at the min pixel 0048 % seg.info(:,6) : second d of the phase parallel to the seg at the min pixel 0049 % seg.info(:,7) and seg_info(:,8) : min and max area of neighboring regions 0050 % seg.info(:,9) and seg_info(:,10) : min and max lengths of the minor axis of the neighboring regions 0051 % seg.info(:,11) and seg_info(:,12) : min and max lengths of the major axis of the neighboring regions 0052 % seg.info(:,11) : length of minor axis 0053 % seg.info(:,12) : length of major axis 0054 % seg.info(:,13) : square of length of major axis 0055 % seg.info(:,16) : max length of region projected onto the major axis 0056 % segment 0057 % seg.info(:,17) : min length of region projected onto the major axis 0058 % segment 0059 % seg.info(:,18) : max length of region projected onto the minor axis 0060 % segment 0061 % seg.info(:,19) : min length of region projected onto the minor axis 0062 % segment 0063 % 0064 % The output images are related by 0065 % mask_cell = mask_bg .* (~segs_good) .* (~segs_3n); 0066 % 0067 % Copyright (C) 2016 Wiggins Lab 0068 % Written by Stella Stylianidou & Paul Wiggins. 0069 % University of Washington, 2016 0070 % This file is part of SuperSegger. 0071 % 0072 % SuperSegger is free software: you can redistribute it and/or modify 0073 % it under the terms of the GNU General Public License as published by 0074 % the Free Software Foundation, either version 3 of the License, or 0075 % (at your option) any later version. 0076 % 0077 % SuperSegger is distributed in the hope that it will be useful, 0078 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0079 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0080 % GNU General Public License for more details. 0081 % 0082 % You should have received a copy of the GNU General Public License 0083 % along with SuperSegger. If not, see <http://www.gnu.org/licenses/>. 0084 0085 % Load the constants from the package settings file 0086 0087 MIN_BG_AREA = CONST.superSeggerOpti.MIN_BG_AREA; 0088 MAGIC_RADIUS = CONST.superSeggerOpti.MAGIC_RADIUS; 0089 MAGIC_THRESHOLD = CONST.superSeggerOpti.MAGIC_THRESHOLD; 0090 CUT_INT = CONST.superSeggerOpti.CUT_INT; 0091 SMOOTH_WIDTH = CONST.superSeggerOpti.SMOOTH_WIDTH; 0092 MAX_WIDTH = CONST.superSeggerOpti.MAX_WIDTH; 0093 A = CONST.superSeggerOpti.A; 0094 verbose = CONST.parallel.verbose; 0095 0096 0097 if ~isfield( CONST.seg, 'segScoreInfo' ) 0098 CONST.seg.segScoreInfo = @segInfoL2; 0099 end 0100 0101 if ~exist('header','var') 0102 header = []; 0103 end 0104 0105 if ~exist('crop_box','var') 0106 crop_box = []; 0107 end 0108 0109 if ~exist('adapt_flag','var') 0110 adapt_flag = 1; 0111 end 0112 0113 if ~exist('mask','var') 0114 mask = []; 0115 end 0116 0117 mask_colonies = []; 0118 0119 % If data structure is passed, rescue the phase image and mask 0120 if isstruct( phaseOrData ) 0121 data = phaseOrData; 0122 phaseOrig = data.phase; 0123 dataFlag = true; 0124 segs_old = data.segs; 0125 0126 % Copy the mask from data into mask if it is empty 0127 if isempty( mask ) 0128 mask = data.mask_bg; 0129 end 0130 else 0131 phaseOrig = phaseOrData; 0132 dataFlag = false; 0133 data = []; 0134 data.phase = phaseOrig; 0135 end 0136 0137 0138 0139 0140 % Initial image smoothing 0141 %this step is necessary to reduce the camera and read noise in the raw 0142 %phase image. Without it, the watershed algorithm will over-segment the 0143 %image. 0144 if all(ismember('100X',CONST.ResFlag)) 0145 phaseNorm = imfilter(phaseOrig,fspecial('disk',1),'replicate'); 0146 else 0147 phaseNorm = phaseOrig; 0148 end 0149 0150 0151 % fix the range, set the max and min value of the phase image 0152 mult_max = 2.5; 0153 mult_min = 0.3; 0154 mean_phase = mean(phaseNorm(:)); 0155 phaseNorm(phaseNorm > (mult_max*mean_phase)) = mult_max*mean_phase; 0156 phaseNorm(phaseNorm < (mult_min*mean_phase)) = mult_min*mean_phase; 0157 0158 0159 % if the size of the meditatrix is even, we get a half pixel shift in the 0160 % position of the mask which turns out to be a probablem later. 0161 f = fspecial('gaussian', 11, SMOOTH_WIDTH); 0162 phaseNormFilt = imfilter(phaseNorm, f,'replicate'); 0163 0164 0165 % Minimum constrast filter to enhance inter-cellular image contrast 0166 [phaseNormFilt,imin,imax] = ag(phaseNormFilt); 0167 magicPhase = magicContrast(phaseNormFilt, MAGIC_RADIUS); 0168 0169 phaseNormUnfilt = (double(phaseOrig)-imin)/(imax-imin); 0170 0171 % C2phase is the Principal curvature 2 of the image without negative values 0172 % it also enhances subcellular contrast. We subtract the magic threshold 0173 % to remove the variation in intesnity within a cell region. 0174 [~,~,~,C2phase] = curveFilter (double(phaseNormFilt),1); 0175 C2phaseThresh = double(uint16(C2phase-MAGIC_THRESHOLD)); 0176 0177 0178 % creates initial background mask by globally thresholding the band-pass 0179 % filtered phase image. We determine the thresholds empirically. 0180 % We use one threshold to remove the background, and another to remove 0181 % the smaller background regions between cells. 0182 if isempty(mask) 0183 % no background making mask 0184 filt_3 = fspecial('gaussian',25, 15); 0185 filt_4 = fspecial('gaussian',5, 1/2); 0186 mask_colonies = makeBgMask(phaseNormFilt, filt_3, filt_4, MIN_BG_AREA, CONST, crop_box); 0187 0188 [~,~,~,~,~,K,~,~] = curveFilter(phaseNormUnfilt, 3 ); 0189 aK = abs(K); 0190 0191 if CONST.superSeggerOpti.remove_debris 0192 mask_colonies = removeDebris(mask_colonies, phaseNormUnfilt, aK, CONST); 0193 end 0194 0195 % remove bright halos from the mask 0196 mask_halos = (magicPhase>CUT_INT); 0197 mask_bg = logical((mask_colonies-mask_halos)>0); 0198 0199 % removes micro-colonies with background level outline intensity - not bright enough 0200 0201 if CONST.superSeggerOpti.remove_microcolonies 0202 mask_bg = intRemoveFalseMicroCol(mask_bg, phaseOrig, CONST); 0203 end 0204 0205 else 0206 mask_bg = mask; 0207 end 0208 0209 if nargin < 3 || isempty(disp_flag) 0210 disp_flag=1; 0211 end 0212 0213 if nargin < 5 || isempty(adapt_flag) 0214 adapt_flag=1; 0215 end 0216 0217 %% Split up the micro colonies into watershed regions to assemble cells 0218 % if data exists, reconstruct ws 0219 if dataFlag 0220 ws = logical(~data.mask_bg+data.segs.segs_3n+data.segs.segs_bad+data.segs.segs_good); 0221 else 0222 % watershed just the cell mask to identify segments 0223 phaseMask = uint8(agd(C2phaseThresh) + 255*(1-(mask_bg))); 0224 ws = 1-(1-double(~watershed(phaseMask,8))).*mask_bg; 0225 0226 0227 if adapt_flag 0228 % If the adapt_flag is set to true (on by default) it watersheds the C2phase 0229 % without using the thershold to identify more segments. It atempts to 0230 % breaks regions that are too big to be cells. This function slows the 0231 % code down, AND slows down the regionOpti code. 0232 0233 wsc = 1- ws; 0234 regs_label = bwlabel( wsc ); 0235 props = regionprops( regs_label, 'BoundingBox','Orientation','MajorAxisLength','MinorAxisLength'); 0236 L2 = [props.MinorAxisLength]; 0237 wide_regions = find(L2 > MAX_WIDTH); 0238 0239 for ii = wide_regions 0240 [xx,yy] = getBB( props(ii).BoundingBox ); 0241 mask_reg = (regs_label(yy,xx)==ii); 0242 0243 c2PhaseReg = double(C2phase(yy,xx)).*mask_reg; 0244 invC2PhaseReg = 1-mask_reg; 0245 ppp = c2PhaseReg+max(c2PhaseReg(:))*invC2PhaseReg; 0246 wsl = double(watershed(ppp)>0); 0247 wsl = (1-wsl).*mask_reg; 0248 0249 % prune added segs by adding just enough to fix the cell width problem 0250 wsl_cc = compConn( wsl, 4 ); 0251 wsl_3n = double(wsl_cc>2); 0252 wsl_segs = wsl-wsl_3n; 0253 wsl_label = bwlabel(wsl_segs,4); 0254 num_wsl_label = max(wsl_label(:)); 0255 wsl_mins = zeros(1,num_wsl_label); 0256 0257 debug_flag = 0; 0258 if debug_flag 0259 backer = 0.5*ag(ppp); 0260 imshow(cat(3,backer,backer,backer + ag(wsl_segs)),[]); 0261 keyboard; 0262 end 0263 0264 for ff = 1:num_wsl_label 0265 wsl_mins(ff) = min(c2PhaseReg(ff==wsl_label)); 0266 end 0267 [wsl_mins, sort_ord] = sort(wsl_mins,'descend'); 0268 0269 wsl_segs_good = wsl_3n; 0270 0271 for ff = sort_ord; 0272 wsl_segs_good = wsl_segs_good + double(wsl_label==ff); 0273 mask_reg_tmp = mask_reg-wsl_segs_good; 0274 if maxMinAxis(mask_reg_tmp) < MAX_WIDTH 0275 break 0276 end 0277 end 0278 ws(yy,xx) = double(0<(ws(yy,xx) + wsl_segs_good)); 0279 end 0280 end 0281 end 0282 0283 %% Remake the data structure 0284 0285 % Determine the "good" and "bad" segments 0286 data = defineGoodSegs(data,ws,phaseNormUnfilt,C2phaseThresh,mask_bg, A,CONST, ~dataFlag ); 0287 data.mask_colonies = mask_colonies; 0288 % copy the existing score into the data structure 0289 if dataFlag 0290 0291 data.segs.score = nan( [size(data.segs.info,1),1] ); 0292 % map the regions 0293 props_tmp = regionprops( data.segs.segs_label, segs_old.segs_label, {'MinIntensity','MaxIntensity'} ); 0294 tmp_min = [props_tmp.MinIntensity]; 0295 tmp_max = [props_tmp.MaxIntensity]; 0296 0297 [y,x] = hist( tmp_max, 0:max(tmp_max)); 0298 good_club = x(y==1); 0299 flagger = and( ismember( tmp_max,good_club), and(tmp_min==tmp_max,tmp_min>0)); 0300 0301 data.segs.score( flagger ) = segs_old.score( tmp_max(flagger) ); 0302 data.segs.score( ~flagger ) = nan; 0303 0304 data.segs.segs_good = segs_old.segs_good; 0305 data.segs.segs_bad = segs_old.segs_bad; 0306 data.segs.segs_3n = segs_old.segs_3n; 0307 0308 %disp( ['New segments lost: ',num2str(sum( ~flagger))] ); 0309 %disp( ['Old segments lost: ',num2str( sum(~ismember(unique(segs_old.segs_label(:)),tmp_max(flagger))))] ); 0310 0311 end 0312 0313 % Calculate and return the final cell mask 0314 data.mask_cell = double((mask_bg - data.segs.segs_good - data.segs.segs_3n)>0); 0315 0316 0317 if disp_flag 0318 figure(1) 0319 clf; 0320 showSegDataPhase(data); 0321 drawnow; 0322 end 0323 0324 0325 0326 0327 end 0328 0329 function Lmax = maxMinAxis(mask) 0330 % maxMinAxis : calculates maximum minor axis length of the regions in the mask. 0331 mask_label = bwlabel(mask); 0332 props = regionprops( mask_label, 'Orientation', 'MinorAxisLength' ); 0333 Lmax = max([props.MinorAxisLength]); 0334 end 0335