multiAssignmentFastOnlyOverlap

 multiAssignmentPairs : links regions in data_c to regions in data_f.
 Each row is assigned to one column only - starting by the minimum
 possible cost and continuing to the next minimum possible cost.

 INPUT :
    (data_c, data_f,CONST, forward, debug_flag)

 OUTPUT :
   [assignments,errorR,totCost,allC,allF]

 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 [assignments,errorR,totCost,allC,allF,dA,revAssign]  = multiAssignmentFastOnlyOverlap (data_c, data_f, CONST, forward, debug_flag)
0002 % multiAssignmentPairs : links regions in data_c to regions in data_f.
0003 % Each row is assigned to one column only - starting by the minimum
0004 % possible cost and continuing to the next minimum possible cost.
0005 %
0006 % INPUT :
0007 %    (data_c, data_f,CONST, forward, debug_flag)
0008 %
0009 % OUTPUT :
0010 %   [assignments,errorR,totCost,allC,allF]
0011 %
0012 % Copyright (C) 2016 Wiggins Lab
0013 % Written by Stella Stylianidou
0014 % University of Washington, 2016
0015 % This file is part of SuperSegger.
0016 %
0017 % SuperSegger is free software: you can redistribute it and/or modify
0018 % it under the terms of the GNU General Public License as published by
0019 % the Free Software Foundation, either version 3 of the License, or
0020 % (at your option) any later version.
0021 %
0022 % SuperSegger is distributed in the hope that it will be useful,
0023 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0024 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0025 % GNU General Public License for more details.
0026 %
0027 % You should have received a copy of the GNU General Public License
0028 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0029 
0030 
0031 DA_MIN = CONST.trackOpti.DA_MIN;
0032 DA_MAX =  CONST.trackOpti.DA_MAX;
0033 OVERLAP_LIMIT_MIN = CONST.trackOpti.OVERLAP_LIMIT_MIN;
0034 
0035 
0036 if ~exist('debug_flag','var') || isempty(debug_flag)
0037     debug_flag = 0;
0038 end
0039 
0040 if ~exist('forward','var') || isempty(forward)
0041     forward = 1;
0042 end
0043 
0044 
0045 if forward
0046     sign = 1;
0047 else
0048     sign = -1;
0049 end
0050 
0051 maxDA = max(sign * DA_MIN,sign * DA_MAX);
0052 minDA = min(sign * DA_MIN,sign * DA_MAX);
0053 
0054 revAssign = [];
0055 assignments = [];
0056 errorR = [];
0057 totCost=[];
0058 allC=[];
0059 allF=[];
0060 
0061 
0062 if ~isempty(data_c)
0063     
0064     if ~isfield(data_c,'regs')
0065         data_c = updateRegionFields (data_c,CONST);
0066     end
0067     
0068     numRegs1 = data_c.regs.num_regs;
0069     assignments = cell( 1, numRegs1);
0070     
0071     errorR = zeros(1,numRegs1);
0072     dA = zeros(1,numRegs1);
0073     
0074     if ~isempty(data_f)
0075         
0076         if ~isfield(data_f,'regs')
0077             data_f = updateRegionFields (data_f,CONST);
0078         end
0079         
0080         areaMat = [data_f.regs.props(:).Area];
0081         numRegs2 = data_f.regs.num_regs;
0082         allC = 1:data_c.regs.num_regs;
0083         allF = 1:data_f.regs.num_regs;
0084         revAssign = cell(1, numRegs2);
0085         totCost = NaN * zeros (numRegs1,numRegs2);
0086         ss = size(data_c.phase);
0087         for ii = 1:numRegs1 % loop through the regions
0088             BB = data_c.regs.props(ii).BoundingBox;
0089             [BBxx,BByy] = getBBpad( BB,ss,5);
0090             tmpc_labels = data_c.regs.regs_label(BByy,BBxx);
0091             maskC =  (tmpc_labels == ii);
0092             areaC =  data_c.regs.props(ii).Area;
0093             
0094             %[maskC,areaC,~] = regProperties (data_c,ii);
0095             
0096             % overlapping regions in f
0097             tmpf_labels = data_f.regs.regs_label(BByy,BBxx);
0098             tmpregs2 = tmpf_labels(maskC);
0099             possibleMapInd = unique(tmpregs2);
0100             possibleMapInd = possibleMapInd(possibleMapInd~=0)'; % remove 0
0101             overlapAreaNorm = zeros(1,numel(possibleMapInd));
0102             
0103             for uu = 1: numel(possibleMapInd)
0104                 % one to one mapping
0105                 regInF = possibleMapInd(uu);
0106                 maskF =  (tmpf_labels == regInF);
0107                 areaF =  data_f.regs.props(regInF).Area;
0108                 overlapMask = maskF(maskC);
0109                 areaOverlap = sum(overlapMask(:)); % area of overlap between jj and ii
0110                 overlapAreaNorm(uu) = areaOverlap/max(areaC,areaF);
0111                 totCost (ii, regInF ) = 1/overlapAreaNorm(uu);
0112             end
0113             
0114             [overlapAreaNormSort, indexesSort] = sort(overlapAreaNorm,'descend');
0115             indexesSort = indexesSort(overlapAreaNormSort>OVERLAP_LIMIT_MIN);
0116             tmpAssign = possibleMapInd(indexesSort);
0117             
0118             % one - to - one mapping
0119             if ~isempty(tmpAssign)
0120                 tmpAreaF = areaMat(tmpAssign);
0121                 totAreaF = cumsum(tmpAreaF);
0122                 areaC = data_c.regs.props(ii).Area;
0123                 tmpdA = (totAreaF-areaC)./totAreaF;
0124                 tmpError = setError(tmpdA);
0125                 
0126                 
0127                 if numel(tmpError) > 0 && tmpError(1) == 0
0128                     assignments{ii} = tmpAssign(1);
0129                     revAssign{tmpAssign(1)} = [revAssign{tmpAssign(1)},ii];
0130                 elseif numel(tmpError) > 1 && tmpError(2) == 0
0131                     assignments{ii} = tmpAssign(1:2);
0132                     revAssign{tmpAssign(1)} = [revAssign{tmpAssign(1)},ii];
0133                     revAssign{tmpAssign(2)} = [revAssign{tmpAssign(2)},ii];
0134                 else
0135                     assignments{ii} = tmpAssign(1);
0136                     revAssign{tmpAssign(1)} = [revAssign{tmpAssign(1)},ii];
0137                 end
0138             end
0139         end
0140         
0141         
0142         [~,minIndxF] = min(totCost,[],2);
0143         [~,minIndxC] = min(totCost,[],1);
0144         cArea = [data_c.regs.props.Area];
0145         fArea = [data_f.regs.props.Area];
0146         
0147         [assignments,revAssign] =fixProblems(assignments,revAssign,minIndxC, cArea, fArea);
0148         [revAssign,assignments] =fixProblems(revAssign,assignments,minIndxF, fArea, cArea);
0149         
0150         dA = changeInArea(assignments, cArea,fArea);
0151         errorR = setError(dA);
0152     end
0153     
0154     % if something assigned to two...
0155     
0156     
0157     
0158     if debug_flag
0159         visualizeLinking(data_c,data_f,assignments);
0160     end 
0161 end
0162 
0163 
0164 
0165     function [assignments, revAssign] =fixProblems (assignments, revAssign, minIndxC, cArea, fArea)
0166         
0167         leftInF = find(cellfun('isempty',revAssign));
0168         
0169         for jj = leftInF
0170             bestAssgnC = minIndxC(jj);
0171             FAlready = assignments{bestAssgnC};
0172             if isempty(FAlready)
0173                 assignments{bestAssgnC} = jj;
0174             else
0175                 revToAlreadyF = revAssign{FAlready};
0176                 areaC = sum(cArea(revToAlreadyF));
0177                 areaFBefore = sum(fArea(FAlready));
0178                 dABefore = (areaFBefore - areaC)/max(areaFBefore,areaC);
0179                 
0180                 if numel(revToAlreadyF) == 2 && ...
0181                         setError(dABefore)>0
0182                     % two assigned to other f - steal one
0183                     areaFjj = fArea(jj);
0184                     newRevToAlreadyF = revToAlreadyF(revToAlreadyF~=bestAssgnC);
0185                     newAreaC = cArea(newRevToAlreadyF);
0186                     areaC = cArea(bestAssgnC);
0187                     newdAjj = (areaFjj - areaC)/max(areaFjj,areaC);
0188                     newdAalreadyF = (areaFBefore - newAreaC)/max(areaFBefore,areaC);;
0189                     if  ~setError(newdAjj) && ...
0190                             ~setError(newdAalreadyF)
0191                         assignments{bestAssgnC} = jj;
0192                         revAssign{jj} = bestAssgnC;
0193                         revAssign{FAlready} = newRevToAlreadyF;
0194                     end
0195                 else
0196                     % see if assigning both to bestAssgnC solves the problem
0197                     tempAssgn = [FAlready,jj];
0198                     areaF = areaFBefore + fArea(jj);
0199                     dAtmp = (areaF - areaC)/max(areaF,areaC);
0200                     if  setError(dABefore) > 0 && ...
0201                             ~setError(dAtmp)
0202                         assignments{bestAssgnC} = tempAssgn;
0203                         revAssign{jj} = bestAssgnC;
0204                     end
0205                 end
0206             end
0207             
0208         end
0209         
0210     end
0211 
0212     function errorR = setError(DA)
0213         errorR = zeros(1, numel(DA));
0214         errorR (DA < minDA ) = 2;
0215         errorR (DA > maxDA) = 3;
0216     end
0217 
0218 
0219 end
0220 
0221 
0222 function dA = changeInArea(assignments, cArea,fArea)
0223 % change in area set for data_c
0224 numRegs1 = size(assignments,2);
0225 dA = nan*zeros(1,numRegs1);
0226 for ll = 1 : numRegs1
0227     tmpAssgn =  assignments{ll};
0228     carea_tmp =  (cArea(ll));
0229     farea_tmp = sum(fArea(tmpAssgn));
0230     dA(ll) = (farea_tmp - carea_tmp) / max(carea_tmp,farea_tmp);
0231 end
0232 
0233 end
0234 
0235 
0236 function visualizeLinking(data_c,data_f,assignments)
0237 figure(1)
0238 clf;
0239 subplot(1,2,1)
0240 imshow(data_c.mask_cell);
0241 subplot(1,2,2)
0242 imshow(data_f.mask_cell);
0243 num_ass = numel(assignments);
0244 randcolor = hsv(256);
0245 markers = {'o','s','d','>','<','^','v','p','h'};
0246 
0247 for c = 1 : num_ass
0248     assF = assignments {c};
0249     if ~isempty(assF)
0250         randomMarker = markers{randi(numel(markers),1)};
0251         randNum = randi(256,1);
0252         color = randcolor(randNum,:);
0253         figure(1);
0254         subplot(1,2,1);
0255         hold on;
0256         plot(data_c.regs.props(c).Centroid(1),data_c.regs.props(c).Centroid(2),[randomMarker,'k'],'MarkerFaceColor',color,'MarkerSize',8);
0257         subplot(1,2,2);
0258         for i = 1 : numel(assF)
0259             hold on;
0260             plot(data_f.regs.props(assF(i)).Centroid(1),data_f.regs.props(assF(i)).Centroid(2),[randomMarker,'k'],'MarkerFaceColor',color,'MarkerSize',8);
0261         end
0262     end
0263 end
0264 end
0265 
0266 function [comboMask,comboArea,comboCentroid] = regProperties (data_c,regNums)
0267 comboCentroid = 0;
0268 comboMask = 0 * (data_c.regs.regs_label);
0269 comboArea = 0;
0270 regNums = regNums(~isnan(regNums));
0271 for ii = 1: numel(regNums)
0272     reg = regNums(ii);
0273     comboCentroid = comboCentroid + data_c.regs.props(reg).Centroid;
0274     comboMask =  comboMask + (data_c.regs.regs_label==reg);
0275     comboArea = comboArea + data_c.regs.props(reg).Area;
0276 end
0277 
0278 comboCentroid = comboCentroid/numel(regNums); % mean centroid
0279 comboMask = (comboMask>0);
0280 end