missingSeg2to1

 missingSeg2to1 : finds missing segment in regC.
 Segments in regC are used that are close to the segment
 between the two regions regR(1) and regR(2) in data_r.
 if a segment is found that fits the requirements data_new is made with
 the new cell_mask and success is returned as true. Else the same data is
 returned and success is false.

 INPUT :
      data_c : current data (seg/err) file.
      regC : numbers of region in current data that needs to be separated
      data_r : reverse data (seg/err) file.
      regR : numbers of regions in reverse data file
      CONST : segmentation parameters

 OUTPUT :
      data_new : data_c with new segment in good_segs and modified cell mask
      success : true if segment was found succesfully.

 Copyright (C) 2016 Wiggins Lab
 Written by Stella Stylianidou
 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_new,success] = missingSeg2to1 (data_c,regC,data_r,regR,CONST)
0002 % missingSeg2to1 : finds missing segment in regC.
0003 % Segments in regC are used that are close to the segment
0004 % between the two regions regR(1) and regR(2) in data_r.
0005 % if a segment is found that fits the requirements data_new is made with
0006 % the new cell_mask and success is returned as true. Else the same data is
0007 % returned and success is false.
0008 %
0009 % INPUT :
0010 %      data_c : current data (seg/err) file.
0011 %      regC : numbers of region in current data that needs to be separated
0012 %      data_r : reverse data (seg/err) file.
0013 %      regR : numbers of regions in reverse data file
0014 %      CONST : segmentation parameters
0015 %
0016 % OUTPUT :
0017 %      data_new : data_c with new segment in good_segs and modified cell mask
0018 %      success : true if segment was found succesfully.
0019 %
0020 % Copyright (C) 2016 Wiggins Lab
0021 % Written by Stella Stylianidou
0022 % University of Washington, 2016
0023 % This file is part of SuperSegger.
0024 %
0025 % SuperSegger is free software: you can redistribute it and/or modify
0026 % it under the terms of the GNU General Public License as published by
0027 % the Free Software Foundation, either version 3 of the License, or
0028 % (at your option) any later version.
0029 %
0030 % SuperSegger is distributed in the hope that it will be useful,
0031 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0033 % GNU General Public License for more details.
0034 %
0035 % You should have received a copy of the GNU General Public License
0036 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0037 
0038 success = false;
0039 data_new = data_c;
0040 debug_flag = 0;
0041 
0042 % need some checks to see if this should happen
0043 longAxis = data_c.regs.info(regC,1);
0044 shortAxis = data_c.regs.info(regC,2);
0045 [xx,yy] = getBBpad(data_c.regs.props(regC).BoundingBox,size(data_c.phase),4);
0046 mask = (data_c.regs.regs_label(yy,xx) == regC);
0047 mask = (mask - data_c.segs.segs_3n(yy,xx))>0;
0048 
0049 % turn on all segments pick best one that divides the area ~ equally?
0050 segsLabel = data_c.segs.segs_label(yy,xx).*mask;
0051 segs_list = unique(segsLabel);
0052 segs_list = segs_list(segs_list~=0);
0053 areaR1 = data_r.regs.props(regR(1)).Area;
0054 areaR2 = data_r.regs.props(regR(2)).Area;
0055 
0056 maskR1 = (data_r.regs.regs_label == regR(1));
0057 maskR2 = (data_r.regs.regs_label == regR(2));
0058 comboMaskR = maskR1+maskR2;
0059 
0060 comboMaskR = comboMaskR(yy,xx);
0061 % find segment between them in data_r
0062 comboMaskRdil = imdilate(comboMaskR, strel('square',3));
0063 comboMaskRerod = imerode(comboMaskRdil, strel('square',3));
0064 separatingSegment = ~comboMaskR.*comboMaskRerod;
0065 
0066 dist = bwdist(separatingSegment);
0067 
0068 if debug_flag
0069     figure(2);
0070     clf;
0071     imshow(cat(3,0.5*ag(data_c.phase) + ag(data_c.regs.regs_label==regC),...
0072         ag(data_r.regs.regs_label==regR(1)),ag(data_r.regs.regs_label==regR(2))));
0073 end
0074 
0075 % do not divide a cell that is too small to be divided.
0076 if numel(segs_list) == 0  || ...
0077         (longAxis < CONST.regionOpti.MIN_LENGTH * 1.5 && shortAxis < .5*CONST.superSeggerOpti.MAX_WIDTH)
0078     return
0079 end
0080 if debug_flag
0081     imshow(segsLabel);
0082 end
0083 
0084 % keep only segments of interest
0085 [minIndex,minRegEScore, minDA,segs_close] = findBestSegs (segsLabel,segs_list,dist,mask,CONST,areaR1,areaR2,data_c.segs.scoreRaw);
0086 
0087 if  ~isempty(minIndex)% && any (minRegEScore) > 0 && minDA < 1.5*CONST.trackOpti.DA_MAX
0088     % a good solution
0089     num_segs = numel(segs_close);
0090     vect = makeVector(minIndex-1,num_segs);
0091     segsAdded = data_c.segs.segs_label * 0;
0092     segsRemoved = data_c.segs.segs_label * 0;
0093     
0094     for kk = 1:num_segs
0095         if  vect(kk)
0096             segsAdded = segsAdded + (segs_close(kk)==data_c.segs.segs_label);
0097         else
0098             segsRemoved = segsRemoved + (segs_close(kk)==data_c.segs.segs_label);
0099         end
0100     end
0101     
0102     segsAdded(yy,xx) = segsAdded(yy,xx) + data_c.segs.segs_3n(yy,xx);
0103     segsAdded(yy,xx) = segsAdded(yy,xx)>0;
0104     data_new.segs.segs_good = data_c.segs.segs_good+segsAdded;
0105     data_new.segs.segs_good = data_c.segs.segs_good-segsRemoved;
0106     data_new.segs.segs_bad = data_c.segs.segs_bad+segsRemoved;
0107     data_new.segs.segs_bad = data_c.segs.segs_bad-segsAdded;
0108     data_new.segs.segs_bad = (data_new.segs.segs_bad>0);
0109     data_new.segs.segs_good = (data_new.segs.segs_good>0);
0110     
0111     data_new.regs.regs_label(segsAdded>0) = 0; % removes segment from regs_label
0112     data_new.mask_cell(segsAdded>0) = 0;
0113     if debug_flag
0114         imshow(cat(3,ag(data_new.mask_cell),ag(data_c.mask_cell),ag(data_c.mask_cell)))
0115     end
0116     success = true;
0117     return;
0118 end
0119 
0120 % resegment only that part of the image
0121 phase = data_c.phase(yy,xx) ;
0122 tmp = superSeggerOpti(phase, [], 1, CONST, 1, '', []);
0123 maskDil =  imdilate(mask, strel('square',2));
0124 newmask = tmp.mask_bg;
0125 newmask(~maskDil) = 0;
0126 if debug_flag
0127     imshow(newmask-tmp.segs.segs_3n-tmp.segs.segs_good-tmp.segs.segs_bad);
0128 end
0129 newmask = (newmask - tmp.segs.segs_3n)>0;
0130 
0131 TmpSegsLabel = tmp.segs.segs_label.*newmask;
0132 segs_list = unique(TmpSegsLabel);
0133 segs_list = segs_list(segs_list~=0);
0134 
0135 [minIndex,minRegEScore, minDA,segs_close] = findBestSegs (TmpSegsLabel,segs_list,dist,newmask,CONST,areaR1,areaR2,tmp.segs.scoreRaw);
0136 
0137 % check if a good solution was found
0138 if  ~isempty(minIndex) %&& any (minRegEScore) > 0 && minDA < 1.5*CONST.trackOpti.DA_MAX
0139     % a good solution
0140     num_segs = numel(segs_close);
0141     vect = makeVector(minIndex-1,num_segs);
0142     segsAdded = TmpSegsLabel * 0;
0143     segsRemoved = TmpSegsLabel * 0;
0144     
0145     for kk = 1:num_segs
0146         if  vect(kk)
0147             segsAdded = segsAdded + (segs_close(kk)==TmpSegsLabel);
0148         else
0149             segsRemoved = segsRemoved + (segs_close(kk)==TmpSegsLabel);
0150         end
0151     end
0152     
0153     finalMask = newmask - segsAdded + segsRemoved;
0154     data_new = data_c;
0155     mask_cell_partial = data_new.mask_cell(yy,xx);
0156     mask_cell_partial((data_c.regs.regs_label(yy,xx) == regC)) = finalMask((data_c.regs.regs_label(yy,xx) == regC));
0157     data_new.mask_cell(yy,xx) = mask_cell_partial;
0158     success = true;
0159     return;
0160 end
0161 
0162 end
0163 
0164 function [minIndex,minRegEScore, minDA,segs_close] = findBestSegs (segsLabel,segs_list,dist,mask,CONST,areaR1,areaR2,rawScores)
0165 verbose = CONST.parallel.verbose;
0166 minimumScore = inf;
0167 DIST_CUT = 5;
0168 dist_segs_c = inf*segs_list;
0169 nsegs = numel(segs_list);
0170 minIndex = [];
0171 minRegEScore = [];
0172 minDA = [];
0173 
0174 for jj = 1:nsegs
0175     dist_segs_c(jj) = min(dist(segsLabel ==segs_list(jj)));
0176 end
0177 
0178 % sort distances and keep only distances < 5
0179 
0180 segs_close = segs_list(dist_segs_c<DIST_CUT);
0181 
0182 if numel(segs_close) > 8 
0183    [~,ord_segs] = sort(dist_segs_c);
0184    segs_close = segs_list(ord_segs(1:8));
0185 end
0186 
0187 % turn on each segment until the regions become more than 1
0188 num_segs = numel( segs_close );
0189 
0190 if num_segs == 0
0191     num_comb = 0;
0192 else
0193     num_comb = 2^num_segs;
0194 end
0195 
0196 if verbose
0197     disp (['Finding missing segment from ', num2str(num_segs),' segments, ', num2str(num_comb), ' combinations']);
0198 end
0199 
0200 for i = 1  : num_comb
0201     
0202     vect = makeVector(i-1,num_segs); % systematic
0203     cell_mask_mod = mask;
0204     segmentMask = mask*0;
0205     for kk = 1:num_segs
0206         if vect(kk)
0207             segmentMask = segmentMask + (segs_close(kk)==segsLabel);
0208         end
0209     end
0210     
0211     % sometimes the problem is that segments are not touching.. dilate the
0212     % segment mask
0213     segmentMask = imdilate(segmentMask, strel('square',2));
0214     cell_mask_mod = (mask - segmentMask) >0;
0215     regs_label_mod = (bwlabel( cell_mask_mod, 4 ));
0216     num_regs_mod   = max(regs_label_mod(:));
0217     
0218     if num_regs_mod == 2 % only if there are two regions
0219         % loop through the regs and get their scores
0220         kk_range = 1:num_regs_mod;
0221         
0222         cell_mod_props = regionprops( regs_label_mod, 'BoundingBox','Orientation','Centroid','Area');
0223         
0224         reg_E = [];
0225         mask_reg = {};
0226         for kk = kk_range;
0227             [xx_,yy_] = getBB( cell_mod_props(kk).BoundingBox);
0228             mask_reg{kk} = (regs_label_mod==kk);
0229             info = CONST.regionScoreFun.props(mask_reg{kk}(yy_,xx_), cell_mod_props(kk));
0230             reg_E(kk) = CONST.regionScoreFun.fun(info,CONST.regionScoreFun.E);
0231         end
0232         
0233         areaChange1 = abs(cell_mod_props(2).Area  - areaR1)/ areaR1 + abs(cell_mod_props(1).Area - areaR2)/areaR2;
0234         areaChange2 = abs(cell_mod_props(2).Area - areaR2)/areaR2 +abs(cell_mod_props(1).Area - areaR1)/areaR1;
0235         
0236         regionScore(i) = 10 * min(areaChange1,areaChange2) + sum(-reg_E)+...
0237             sum((1-2*vect).*rawScores(segs_close)')*CONST.regionOpti.DE_norm;
0238         
0239         if regionScore(i) < minimumScore
0240             minimumScore =  regionScore(i) ;
0241             minIndex = i ;
0242             minRegEScore = reg_E;
0243             minDA = min(areaChange1,areaChange2)/2; % for two regions..
0244         end
0245     end
0246 end
0247 
0248 
0249 end
0250 
0251 function vect = makeVector( nn, n )
0252 vect = zeros(1,n);
0253 for i=n-1:-1:0;
0254     vect(i+1) = floor(nn/2^i);
0255     nn = nn - vect(i+1)*2^i;
0256 end
0257 end