Description of perRegionOpti

0001 function [data] = perRegionOpti( data, disp_flag, CONST,header)
0002 % perRegionOpti : Segmentaion optimization using region characteristics.
0003 % It attempts to improve the region score by turning off on and off segments
0004 % of regions with bad scores. It uses systematic method, or simulated
0005 % anneal, according to the nudmber of segments to be considered.
0006 %
0007 % INPUT :
0008 %       data : data with segs field (.err data or .trk data)
0009 %       dissp : display flag
0010 %       CONST : segmentation constants
0011 %       header : information string
0012 % OUTPUT :
0013 %       data : data structure with modified segments
0014 %
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 
0035 MIN_LENGTH = CONST.regionOpti.MIN_LENGTH;
0036 CutOffScoreHi = CONST.regionOpti.CutOffScoreHi;
0037 MAX_NUM_RESOLVE  = CONST.regionOpti.MAX_NUM_RESOLVE;
0038 MAX_NUM_SYSTEMATIC = CONST.regionOpti.MAX_NUM_SYSTEMATIC;
0039 NUM_INFO = CONST.regionScoreFun.NUM_INFO;
0040 E = CONST.regionScoreFun.E;
0041 
0042 minGoodRegScore = CONST.regionOpti.minGoodRegScore ;
0043 neighMaxScore = CONST.regionOpti.neighMaxScore;
0044 verbose = CONST.parallel.verbose;
0045 
0046 if ~exist('header','var')
0047     header = [];
0048 end
0049 
0050 if nargin < 2 || isempty(disp_flag);
0051     disp_flag = 1;
0052 end
0053 
0054 debug_flag = 0;
0055 
0056 segsLabelAll = data.segs.segs_label;
0057 segs_3n    = data.segs.segs_3n;
0058 above_Hi_ind = find(data.segs.scoreRaw > CutOffScoreHi);
0059 below_Hi_ind = find(data.segs.scoreRaw <= CutOffScoreHi);
0060 segs_HiGood = double(ismember(segsLabelAll, above_Hi_ind));
0061 segs_good = segs_HiGood;
0062 segs_bad = double(ismember(segsLabelAll, below_Hi_ind));
0063 ss = size(segs_3n);
0064 
0065 % remaining segs allowed to be tweaked
0066 segsLabelMod = data.segs.segs_label;
0067 segsLabelMod(logical(segs_3n+segs_good)) = 0;
0068 
0069 % remake mask with best guessed regions : segs_3n and high segs_good
0070 mask_regs = double((data.mask_bg-segs_3n-segs_good)>0);
0071 data.regs.regs_label = (bwlabel( mask_regs, 4 ));
0072 data.regs.props = regionprops( data.regs.regs_label, ...
0073     'BoundingBox','Orientation','Centroid','Area');
0074 data.regs.num_regs = max( data.regs.regs_label(:) );
0075 data.regs.score  = ones( data.regs.num_regs, 1 );
0076 data.regs.scoreRaw = ones( data.regs.num_regs, 1 );
0077 data.regs.info = zeros( data.regs.num_regs, NUM_INFO );
0078 
0079 for ii = 1:data.regs.num_regs
0080     [xx,yy] = getBBpad( data.regs.props(ii).BoundingBox, ss, 1);
0081     mask = data.regs.regs_label(yy,xx)==ii;
0082     data.regs.info(ii,:) = CONST.regionScoreFun.props( mask, data.regs.props(ii) );
0083     data.regs.scoreRaw(ii) = CONST.regionScoreFun.fun(data.regs.info(ii,:), E);
0084     data.regs.score(ii) = data.regs.scoreRaw(ii) > 0;
0085 end
0086 
0087 if verbose
0088     disp([header, 'rO: Got ',num2str(data.regs.num_regs),' regions.']);
0089 end
0090 
0091 % get regions with bad scores
0092 regs_label = data.regs.regs_label;
0093 badReg = find(data.regs.scoreRaw < minGoodRegScore);
0094 props = data.regs.props;
0095 
0096 numBadRegions = size(badReg,1);
0097 if verbose
0098     disp([header, 'rO: Possible segments to be tweaked : ',num2str(numel(unique(segsLabelMod))-1),'.']);
0099     disp([header, 'rO: Optimizing ',num2str(numBadRegions),' regions.']);
0100 end
0101 
0102 % list of already tweaked segments
0103 goodSegList = [];
0104 badSegList = [];
0105 
0106 
0107 while ~isempty(badReg)
0108     ii = badReg(1);
0109     
0110     % get padded box
0111     originalBBbox = data.regs.props(ii).BoundingBox;
0112     [xx_,yy_] = getBBpad( originalBBbox, ss, 5);
0113     cellMask = (regs_label(yy_,xx_) == ii);
0114     
0115     % get neigbors with scores below neighMaxScore
0116     neighborIds = unique(regs_label(yy_,xx_));
0117     neighborIds = neighborIds(neighborIds~=0);
0118     neighbordIds = neighborIds(data.regs.scoreRaw(neighborIds) < neighMaxScore);
0119     cellMaskDil = imdilate(cellMask, strel('square',3));
0120     combinedCellMask = regs_label *0;
0121     
0122     for kk = 1 : numel(neighborIds)
0123         combinedCellMask = combinedCellMask + (regs_label == neighborIds(kk));
0124     end
0125     
0126     % add all segments inside the original cell / touching neighbors
0127     if  data.regs.info(ii,1) < MIN_LENGTH % for a small cell add all segments
0128         segs_list = unique(cellMaskDil.*segsLabelAll(yy_,xx_));
0129     else % add only segments to be considered
0130         segs_list = unique(cellMaskDil.*segsLabelMod(yy_,xx_));
0131     end
0132     segs_list = segs_list(segs_list~=0);
0133     
0134     % add segments to the mask
0135     for kk = 1 :numel(segs_list)
0136         combinedCellMask = combinedCellMask + (segsLabelAll==segs_list(kk));
0137     end
0138     combinedCellMask = combinedCellMask>0;
0139     
0140     % remake regions
0141     regCombLabels = bwlabel(combinedCellMask);
0142     combRegProps = regionprops(regCombLabels,'BoundingBox','Orientation','Centroid','Area');
0143     
0144     % find most overlapping region
0145     origCellMask = (regs_label == ii);
0146     %imshow(cat(3,0.5*ag(mask_regs),ag(origCellMask),ag(origCellMask)));
0147     overlapArea = zeros(1,numel(combRegProps));
0148     for j = 1 : numel(combRegProps)
0149         mask2 = (regCombLabels == j);
0150         overlap = origCellMask & mask2;
0151         overlapArea(j) = sum(overlap(:)) / data.regs.props(ii).Area;
0152     end
0153     
0154     
0155     ind = find(overlapArea ==1);
0156     
0157     % bounding box of overlapping region
0158     [xx,yy] = getBBpad(combRegProps(ind).BoundingBox,ss,3);
0159     combMask = (regCombLabels(yy,xx) == ind);
0160     
0161     % get indices of neighbors we are checking and remove them from
0162     % checklist
0163     neighborsToCheck = unique(regs_label.* (regCombLabels == ind));
0164     badReg = setdiff (badReg,neighborsToCheck);
0165     
0166     % dilate and erode mask to add more segments.
0167     combMaskDil = imdilate(combMask, strel('square',2));
0168     combMaskErode = imerode(combMaskDil, strel('square',2));
0169     segs_list_extra = unique(combMaskErode.*segsLabelMod(yy,xx));
0170     segs_list = [segs_list;segs_list_extra];
0171     segs_list = unique(segs_list);
0172     segs_list = segs_list(segs_list~=0);
0173     
0174     goodChecked = segs_list(ismember(segs_list,goodSegList));
0175     badChecked = segs_list(ismember(segs_list,badSegList));
0176     
0177     %remove from mask already checked segments
0178     for kk = 1 : numel(goodChecked)
0179         combMask = combMask - (segsLabelMod(yy,xx)==goodChecked(kk));
0180     end
0181     
0182     
0183     for kk = 1 : numel(badChecked)
0184         combMask = combMask + (segsLabelMod(yy,xx)==badChecked(kk));
0185     end
0186     
0187     combMask = double(combMask>0);
0188     %keep only not checked segs
0189     if ~isempty(goodChecked)
0190         segs_list =  segs_list(~ismember(segs_list,goodChecked));
0191     end
0192     if ~isempty(badChecked)
0193         segs_list =  segs_list(~ismember(segs_list,badChecked));
0194     end
0195     
0196     
0197     num_segs = numel(segs_list);
0198     segmentMask = 0 *combinedCellMask;
0199     
0200     for kk = 1 :num_segs
0201         segmentMask = segmentMask + (segsLabelAll==segs_list(kk));
0202     end
0203     
0204     
0205     
0206     if isempty(segs_list)
0207         [vect] = [];
0208     elseif numel(segs_list) > MAX_NUM_RESOLVE % use raw score
0209         if verbose
0210             disp([header, 'rO: Too many regions to analyze (',num2str(num_segs),').']);
0211         end
0212         [vect] = data.segs.scoreRaw(segs_list)>0;
0213     elseif numel(segs_list) > MAX_NUM_SYSTEMATIC % use simulated anneal
0214         if verbose
0215             disp([header, 'rO: Simulated Anneal : (',num2str(num_segs),' segments).']);
0216         end
0217         [vect,~] = simAnnealMap( segs_list, data, combMask, xx, yy, CONST, 0);
0218     else % use systematic
0219         if verbose
0220             disp([header, 'rO: Systematic : (',num2str(num_segs),' segments).']);
0221         end
0222         [vect,~] = systematic( segs_list, data, combMask, xx, yy, CONST);
0223     end
0224     
0225     if ~isempty(segs_list)
0226         goodSegList = [goodSegList;segs_list(logical(vect))];
0227         badSegList = [badSegList;segs_list(~vect)];
0228     end
0229     
0230     
0231     if debug_flag &&  ~isempty(segs_list)
0232         % shows the final optimized segments for the combined mask
0233         mod = mask_regs*0;
0234         mod(yy,xx) = combMask;
0235         segment_mask = mask_regs*0;
0236         % segment_mask
0237         for kk = 1:num_segs
0238             mod = mod - vect(kk)*(segs_list(kk)==segsLabelAll);
0239             segment_mask = segment_mask + vect(kk)*(segs_list(kk)==segsLabelAll);
0240         end
0241         
0242         
0243         backer = 0.5*ag(mask_regs);
0244         segment_mask_small = ag(segment_mask(yy,xx));
0245         backer_small = (backer(yy,xx));
0246         backer_small = ag(((backer_small) - ag(combMask))>0);
0247         all_segs = ag(ismember(segsLabelAll,segs_list));
0248         all_segs_small = all_segs(yy,xx);
0249         combMask_bef = data.mask_cell(yy,xx);
0250         figure (2);
0251         subplot(1,2,1)
0252         imshow(cat(3,0.5*(backer_small),backer_small + ag(combMask_bef), ...
0253             ag(combMask_bef) + (backer_small)))
0254         subplot(1,2,2)
0255         imshow(cat(3,0.3*(backer_small)+ag(segment_mask_small)...
0256             ,0.2*ag(all_segs_small-segment_mask_small) +  0.3*backer_small + 0.5*ag(combMask>0), ...
0257             0.5* ag(combMask>0) +  0.3*(backer_small)));
0258         keyboard;
0259     end
0260 end
0261 
0262 if ~isempty(goodSegList)
0263     segs_good = (ismember(segsLabelAll,goodSegList)+segs_good)>0;
0264     segs_bad = (segs_bad - ismember(segsLabelAll,goodSegList))>0;
0265     data.segs.score(goodSegList) = 1;
0266 end
0267 
0268 if ~isempty(badSegList)
0269     segs_good = (segs_good - ismember(segsLabelAll,badSegList))>0;
0270     segs_bad = (segs_bad + ismember(segsLabelAll,badSegList ))>0;
0271     data.segs.score(badSegList) = 0;
0272 end
0273 
0274 data.mask_cell = double((data.mask_bg-segs_3n-segs_good)>0);
0275 data.segs.segs_good = segs_good;
0276 data.segs.segs_bad = segs_bad;
0277 
0278 % update region fields using new mask
0279 data = intMakeRegs( data, CONST );
0280 
0281 if disp_flag
0282     figure(1);
0283     cell_mask = data.mask_cell;
0284     back = double(0.7*ag( data.phase ));
0285     outline = imdilate( cell_mask, strel( 'square',3) );
0286     outline = ag(outline-cell_mask);
0287     % imshow does not work for parallel - use image
0288     image(uint8(cat(3,back + 0.7*double(ag(data.segs.segs_good)) + double(ag(data.segs.segs_3n))...
0289         ,back,...
0290         back+ 0.2*double(ag(~cell_mask)-outline))));
0291     drawnow;
0292     axis equal tight;
0293 end
0294 
0295 
0296 end
0297
perRegionOpti

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
