Description of regionOpti

0001 function [data] = regionOpti( data, disp_flag, CONST,header)
0002 % regionOpti : Segmentaion optimization using region characteristics.
0003 % It turns off on and off segments that have scores between two values in the
0004 % constants (CONST.regionOpti.CutOffScoreHi and CONST.regionOpti.CutOffScoreLo)
0005 % And uses systematic method, or simulated anneal, to find the optimal segments
0006 % configuration.
0007 %
0008 % INPUT :
0009 %       data : data with segs field (.err data or .trk data)
0010 %       disp_flag : display flag
0011 %       CONST : segmentation constants
0012 %       header : information string
0013 % OUTPUT :
0014 %       data : data structure with modified segments
0015 %
0016 % Copyright (C) 2016 Wiggins Lab
0017 % Written by Stella Styliandou & Paul Wiggins.
0018 % University of Washington, 2016
0019 % This file is part of SuperSegger.
0020 %
0021 % SuperSegger is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % SuperSegger is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0033 
0034 MAX_WIDTH = CONST.superSeggerOpti.MAX_WIDTH;
0035 MIN_LENGTH = CONST.regionOpti.MIN_LENGTH;
0036 CutOffScoreHi = 30; %CONST.regionOpti.CutOffScoreHi;
0037 CutOffScoreLo = -30; %;CONST.regionOpti.CutOffScoreLo;
0038 MAX_NUM_RESOLVE = CONST.regionOpti.MAX_NUM_RESOLVE;
0039 MAX_NUM_SYSTEMATIC = CONST.regionOpti.MAX_NUM_SYSTEMATIC;
0040 CONST.regionOpti.Emin  = .2;
0041 DE_norm = CONST.regionOpti.DE_norm;
0042 verbose = CONST.parallel.verbose;
0043 
0044 if ~exist('header')
0045     header = [];
0046 end
0047 
0048 if nargin < 2 || isempty('dispp');
0049     disp_flag = 1;
0050 end
0051 
0052 
0053 % Turn on and off segs outside the cutoff.
0054 segs_label = data.segs.segs_label;
0055 segs_3n = data.segs.segs_3n;
0056 segs_bad = 0*data.segs.segs_3n;
0057 segs_good  = segs_bad;
0058 segs_good_off  = segs_bad;
0059 ss = size(segs_3n);
0060 num_segs_ = numel(data.segs.score);
0061 
0062 % segments above the high cutoff are added to segs_3n (hard segments,
0063 % always on) and segments below the low cutoff are added to segs_bad,
0064 % segs label contains only segments that can be switched on / off.
0065 % Both are removed from the local copy of segs_label, so they can not
0066 % be switched on and off.
0067 above_Hi_ind = find(data.segs.scoreRaw > CutOffScoreHi);
0068 below_Lo_ind = find(data.segs.scoreRaw < CutOffScoreLo);
0069 segs_3n = segs_3n + double(ismember(segs_label, above_Hi_ind)); % high in 3n
0070 segs_bad = double(ismember(segs_label, below_Lo_ind)); % low in bad
0071 segs_label(logical(segs_3n+segs_bad)) = 0; % rest remain in segs_label
0072 
0073 mask_regs = double((data.mask_bg-segs_3n)>0);
0074 regs_label = (bwlabel( mask_regs, 4 ));
0075 regs_props = regionprops( regs_label, 'BoundingBox','Orientation' );
0076 num_regs   = max(regs_label(:));
0077 segs_added = [];
0078 if verbose
0079 disp([header, 'rO: Got ',num2str(num_regs),' regions.']);
0080 end
0081 % Find short regions and add surrounding segments to segs_added
0082 for ii = 1:num_regs
0083     
0084     [xx,yy] = getBBpad(regs_props(ii).BoundingBox,ss,2);
0085     tmp_mask = (regs_label(yy,xx)==ii);
0086     % calculates long and short axis of region
0087     [L1,L2] = makeRegSize (tmp_mask, regs_props(ii));
0088     debug_flag = 0;
0089     
0090     if debug_flag
0091         % image of phase and regions (green). Current region is shown in yellow
0092         figure;
0093         clf;
0094         imshow( cat(3,ag(regs_label==ii),ag(regs_label>0),ag(data.phase)), [])
0095         disp([num2str(L1),', ',num2str(L2)]);
0096     end
0097     
0098      
0099     % if region is shorter than MIN_LENGTH it adds the hard segments inside the
0100     % region in segs_added to be switched on / off.
0101     if L1 < MIN_LENGTH;
0102         tmp_mask = imdilate(tmp_mask, strel('square',3));
0103         tmp_added = unique( tmp_mask.*data.segs.segs_label(yy,xx).*segs_3n(yy,xx));
0104         tmp_added = tmp_added(logical(tmp_added));
0105         tmp_added = reshape(tmp_added,1,numel(tmp_added));
0106         segs_added = [segs_added,tmp_added];
0107     end
0108     
0109 end
0110 
0111 % segs_added are surrounding segments of small regions
0112 segs_added = unique(segs_added);
0113 segs_added_ = ismember( data.segs.segs_label, segs_added); % image of segs_added segments
0114 segs_3n(segs_added_) = 0; % removes segs_added_ from 3n
0115 
0116 % adding them to the local copy of segs_label
0117 segs_label(segs_added_) = data.segs.segs_label(segs_added_); 
0118 
0119 
0120 % mask_regs is the super mask of all boundaries + segments that are
0121 % permanently on
0122 mask_regs = double((data.mask_bg-segs_3n)>0);
0123 regs_label = (bwlabel( mask_regs, 4 )); % labels these regions.
0124 regs_props = regionprops( regs_label, 'BoundingBox','Orientation'  );
0125 num_regs   = max( regs_label(:));
0126 
0127 rs_list = cell(1,num_regs);
0128 ss = size(data.phase);
0129 
0130 for ii = 1:num_regs
0131     [xx,yy] = getBBpad(regs_props(ii).BoundingBox,ss,2);
0132     cell_mask = (regs_label(yy,xx) == ii);
0133     
0134     % get the names of the remaining segments that are in this region.
0135     segs_list = unique( cell_mask.*segs_label(yy,xx));
0136     
0137     % get rid of zero and make sure the thing is a row vector.
0138     segs_list = segs_list(logical(segs_list));
0139     segs_list = reshape(segs_list,1,numel(segs_list));   
0140     rs_list{ii} = segs_list;
0141     
0142  
0143     % Turn on segments who would help resolve cells that are too wide
0144     % First turn everything on in the seg_list and check to make sure that
0145     % the regions are small enough.
0146     tmp_segs = cell_mask-ismember( segs_label(yy,xx),segs_list );
0147     tmp_segs = double(tmp_segs>0);
0148     tmp_label = (bwlabel( tmp_segs, 4 ));
0149     tmp_props = regionprops( tmp_label, 'BoundingBox','Orientation'  );
0150     num_tmp   = max( tmp_label(:));    
0151     segs_added = [];
0152     
0153     for ff = 1:num_tmp
0154         
0155         tmp_mask = (tmp_label==ff);
0156         [L1,L2] = makeRegSize( tmp_mask, tmp_props(ff) );
0157         
0158         if L2 > MAX_WIDTH; % break down things larger than max width
0159             tmp_added = unique( tmp_mask.*data.segs.segs_label(yy,xx));
0160             tmp_added = tmp_added(logical(tmp_added));
0161             tmp_added = reshape(tmp_added,1,numel(tmp_added));
0162             segs_added = [segs_added,tmp_added];
0163         end
0164     end
0165     
0166     segs_list = unique([segs_list,segs_added]);
0167 
0168     
0169     if isempty(segs_list)
0170         [vect] = [];
0171     elseif numel(segs_list) > MAX_NUM_RESOLVE % use raw score
0172         if verbose
0173         disp([header, 'rO: Too many regions to analyze (',num2str(numel(segs_list)),').']);
0174         end
0175         [vect] = data.segs.scoreRaw(segs_list)>0;
0176     elseif numel(segs_list) > MAX_NUM_SYSTEMATIC % use simulated anneal
0177          if verbose
0178         disp([header, 'rO: Simulated Anneal : (',num2str(numel(segs_list)),' segments).']);
0179          end
0180         debug_flag = 0;
0181         [vect] = simAnnealMap( segs_list, data, cell_mask, xx, yy, CONST, debug_flag);
0182     else % use systematic
0183         if verbose 
0184         disp([header, 'rO: Systematic : (',num2str(numel(segs_list)),' segments).']);    
0185         end
0186         %tic;
0187         [vect] = systematic( segs_list, data, cell_mask, xx, yy, CONST);
0188        % toc;
0189     end
0190     
0191     num_segs = numel(segs_list);
0192     
0193     try
0194         segs_good(yy,xx) = segs_good(yy,xx) + ismember( data.segs.segs_label(yy,xx), segs_list(logical(vect)));
0195         segs_good_off(yy,xx) = segs_good_off(yy,xx) + ismember( data.segs.segs_label(yy,xx), segs_list(~vect));
0196     catch ME
0197         printError(ME);
0198     end
0199     
0200 end
0201 
0202 data.mask_cell = double((data.mask_bg-segs_3n-segs_good)>0);
0203 
0204 % reset the seg scores incase you want to use this with segsManage
0205 segs_on_ind = unique((~data.mask_cell).*data.segs.segs_label);
0206 segs_on_ind = segs_on_ind(logical(segs_on_ind));
0207 data.segs.score(segs_on_ind) = 1;
0208 segs_off_ind = unique((data.mask_cell).*data.segs.segs_label);
0209 segs_off_ind = segs_off_ind(logical(segs_off_ind));
0210 data.segs.score(segs_off_ind) = 0;
0211 cell_mask = data.mask_cell;
0212 data.segs.segs_good   = double(data.segs.segs_label>0).*double(~data.mask_cell);
0213 data.segs.segs_bad   = double(data.segs.segs_label>0).*data.mask_cell;
0214 
0215 data = intMakeRegs( data, CONST );
0216 
0217 if disp_flag
0218     back = double(0.7*ag( data.phase ));
0219     outline = imdilate( cell_mask, strel( 'square',3) );
0220     outline = ag(outline-cell_mask);
0221     segs_never = ((segs_bad-segs_good_off-segs_good)>0);
0222     segs_tried = ((segs_good_off + segs_good)>0);
0223     
0224     try
0225         figure(1);
0226         clf;
0227         
0228         imshow(uint8(cat(3,back + 1*double(outline),...
0229             back + 0.3*double(ag(segs_tried)),...
0230             back + 0.3*double(ag(segs_tried)).*double(segs_good_off) + 0.2*double(ag(~cell_mask)-outline) + 0.5*double(ag(segs_never)))));
0231     catch ME
0232         printError(ME);
0233     end
0234     drawnow;
0235     
0236 end
0237 end
0238
regionOpti

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
