multiAssignmentSparse

 multiAssignmentSparse : assigns regions in data_c to regions in data_f.
 Uses a combination of area overlap, centroid distance, and outward push
 in colonies. Regions are assigned one-to-one or one-to-pair or
 pair-to-one. Each cell is assigned only once - starting by the
 minimum possible cost and continuing to the next minimum possible cost.
 Uses a sparse matrix to save the costs.

 INPUT :
    data_c : current frame file
    data_f : forward frame file
    CONST : segmentation parameters
    forward : 1 for forward direction (e.g current to forward), 0 for
    reverse
    debug_flag : 1 to display assignment result.

 OUTPUT :
   assignments : cell matrix. cell of region c is assigned id of region f 
   errorR : matrix with 2 (DA<MIN) or 3(DA>MAX) if error, 0 if no error
   totCost : cost matrix
   indexC : region ids in data_c for cost matrix
   indexF : region ids in data_f for cost matrix
   dA :  area change
   revAssign : cell matrix with reverse assignments (from ids in f to c)

 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,indexC,indexF,dA,revAssign]  = multiAssignmentSparse ...
0002     (data_c, data_f,CONST, forward, debug_flag)
0003 % multiAssignmentSparse : assigns regions in data_c to regions in data_f.
0004 % Uses a combination of area overlap, centroid distance, and outward push
0005 % in colonies. Regions are assigned one-to-one or one-to-pair or
0006 % pair-to-one. Each cell is assigned only once - starting by the
0007 % minimum possible cost and continuing to the next minimum possible cost.
0008 % Uses a sparse matrix to save the costs.
0009 %
0010 % INPUT :
0011 %    data_c : current frame file
0012 %    data_f : forward frame file
0013 %    CONST : segmentation parameters
0014 %    forward : 1 for forward direction (e.g current to forward), 0 for
0015 %    reverse
0016 %    debug_flag : 1 to display assignment result.
0017 %
0018 % OUTPUT :
0019 %   assignments : cell matrix. cell of region c is assigned id of region f
0020 %   errorR : matrix with 2 (DA<MIN) or 3(DA>MAX) if error, 0 if no error
0021 %   totCost : cost matrix
0022 %   indexC : region ids in data_c for cost matrix
0023 %   indexF : region ids in data_f for cost matrix
0024 %   dA :  area change
0025 %   revAssign : cell matrix with reverse assignments (from ids in f to c)
0026 %
0027 % Copyright (C) 2016 Wiggins Lab
0028 % Written by Stella Stylianidou
0029 % University of Washington, 2016
0030 % This file is part of SuperSegger.
0031 %
0032 % SuperSegger is free software: you can redistribute it and/or modify
0033 % it under the terms of the GNU General Public License as published by
0034 % the Free Software Foundation, either version 3 of the License, or
0035 % (at your option) any later version.
0036 %
0037 % SuperSegger is distributed in the hope that it will be useful,
0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0040 % GNU General Public License for more details.
0041 %
0042 % You should have received a copy of the GNU General Public License
0043 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0044 
0045 
0046 global str8
0047 str8 = strel('square',5);
0048 DA_MIN = CONST.trackOpti.DA_MIN;
0049 DA_MAX =  CONST.trackOpti.DA_MAX;
0050 
0051 if ~exist('debug_flag','var') || isempty(debug_flag)
0052     debug_flag = 0;
0053 end
0054 
0055 if ~exist('forward','var') || isempty(forward)
0056     forward = 1;
0057 end
0058 
0059 if forward
0060     sign = 1;
0061 else
0062     sign = -1;
0063 end
0064 
0065 maxDA = max(sign * DA_MIN,sign * DA_MAX);
0066 minDA = min(sign * DA_MIN,sign * DA_MAX);
0067 
0068 
0069 
0070 revAssign = [];
0071 assignments = [];
0072 errorR = [];
0073 totCost=[];
0074 indexC=[];
0075 indexF=[];
0076 noOverlap = 0.0001;
0077 centroidWeight = 5;
0078 areaFactor = 20;
0079 areaChangeFactor = 100;
0080 outwardMotFactor = 1/10;
0081 dA = [];
0082 
0083 if ~isempty(data_c)
0084     if ~isfield(data_c,'regs')
0085         data_c = updateRegionFields (data_c,CONST);
0086     end
0087     ss = size(data_c.phase);
0088     
0089     numRegs1 = data_c.regs.num_regs;
0090     assignments  = cell( 1, numRegs1);
0091     errorR = zeros(1,numRegs1);
0092     dA = nan*zeros(1,numRegs1);
0093     
0094     if ~isempty(data_f)
0095         if ~isfield(data_f,'regs')
0096             data_f = updateRegionFields (data_f,CONST);
0097         end
0098         
0099         numRegs2 = data_f.regs.num_regs;
0100         regsInC = 1:data_c.regs.num_regs;
0101         regsInF = 1:data_f.regs.num_regs;
0102         
0103         idsC(1,regsInC) = regsInC;
0104         idsC(2,regsInC) = NaN;
0105         
0106         idsF(1,regsInF) = regsInF;
0107         idsF(2,regsInF) = NaN;
0108         
0109         % colony calculation
0110         maskBgFill= imfill(data_c.mask_bg,'holes');
0111         colony_labels = bwlabel(maskBgFill);
0112         colony_props = regionprops( colony_labels,'Centroid','Area');
0113         
0114         % find possible pairs
0115         [pairsF,~] = findNeighborPairs (data_f, numRegs2, regsInF,[]);
0116         allF = [idsF,pairsF];
0117         
0118         [pairsC,neighF] = findNeighborPairs (data_c, numRegs1, regsInC, data_f);
0119         allC = [idsC,pairsC];
0120         
0121         % initialize
0122         
0123         
0124         indexC =  NaN * ones(2,size(allC,2)*10);
0125         indexF = indexC;
0126         totCost =  NaN * ones(1,size(allC,2)*10);
0127         areaOverlapCost = totCost;
0128         areaChange = totCost;
0129         centroidCost = totCost;
0130         areaOverlapTransCost = totCost;
0131         outwardMot = totCost;
0132         distFromColonyMat = totCost;
0133         goodOneToOne = zeros(1,size(allC,2));
0134         distFromColn = zeros(1,size(allC,2));
0135         counter = 1;
0136         
0137         for ii = 1:size(allC,2) % loop through the regions
0138             % ind : list of regions that overlap with region ii in data 1
0139             
0140          
0141             % if it already has a good mapping don't bother..
0142             cRegs = allC(:,ii);
0143             isSingleRegC = sum(isnan(cRegs)); % has a nan
0144             
0145             % condition for good one to one to mapping
0146             alreadyFoundOneToOne = ~isSingleRegC  && (goodOneToOne(cRegs(1)) ...
0147                 || goodOneToOne(cRegs(2))) ;
0148             
0149             % only check pairs if not good one-to-one mapping
0150             if  ~alreadyFoundOneToOne %&& regionExists
0151                 
0152                 [BB_c_xx,BB_c_yy] = getBoxLimits (data_c,cRegs);
0153                 [maskC,areaC,centroidC] = regProperties (data_c,cRegs,BB_c_xx,BB_c_yy);
0154                 
0155                 % colony it belongs to
0156                 colony_labels_temp = colony_labels(BB_c_yy,BB_c_xx);
0157                 colOverlap = colony_labels_temp(maskC);
0158                 if sum(colOverlap(:)) == 0
0159                     distFromColony = [0 ,0];
0160                     distFromColn (ii) = sqrt(sum(distFromColony.^2));
0161                 else
0162                     colonyId =  max(colOverlap);
0163                     distFromColony = centroidC - colony_props(colonyId).Centroid;
0164                     distFromColn (ii) = sqrt(sum(distFromColony.^2));
0165                 end
0166                 
0167                 % dilate mask and get adjacent regions within dilated area
0168                 nonNanCregs = cRegs(~isnan(cRegs));
0169                 tmpregs2 = [neighF{[nonNanCregs]}];
0170                 possibleMapInd = unique(tmpregs2);
0171                 possibleMapInd = possibleMapInd(possibleMapInd~=0)'; % remove 0
0172                 
0173                 startcounter = counter;
0174                 for yy = 1:numel(possibleMapInd)
0175                     % one to one mapping
0176                     idF = possibleMapInd(yy);
0177                     [maskF,areaF,centroidF] = regProperties (data_f,idF,BB_c_xx,BB_c_yy);
0178                     indexC(:,counter) = cRegs;
0179                     indexF(1,counter) = idF;
0180                     overlapMask = maskF(maskC);
0181                     areaOverlap = sum(overlapMask(:)); % area of overlap between jj and ii
0182                     areaOverlapCost(counter) = areaOverlap/areaC;
0183                     
0184                     if  (areaOverlapCost(counter) == 0)
0185                         areaOverlapCost(counter) = noOverlap;
0186                     end
0187                     
0188                     displacement = centroidC - centroidF;
0189                     centroidCost(counter) = sqrt(sum((displacement).^2));
0190                     if centroidCost(counter) == 0
0191                         outwardMot(counter) = 0;
0192                     else
0193                         outwardMot(counter) = (distFromColony*displacement')/centroidCost(counter);
0194                     end
0195                     
0196                     distFromColonyMat (counter) = exp(-distFromColn(ii)/outwardMotFactor);
0197                     
0198                     areaChange(counter) = (areaF - areaC)/(areaC);
0199                     
0200                     % moved area
0201                     offset = round(displacement);
0202                     maskOut = imshift(maskF,offset);
0203                     maskOut = maskOut(maskC);
0204                     areaOverlapTrns = sum(maskOut(:));
0205                     areaOverlapTransCost(counter) = areaOverlapTrns/areaC;
0206                     
0207                     if (areaOverlapTransCost(counter) == 0)
0208                         areaOverlapTransCost(counter) = noOverlap;
0209                     end
0210                     counter = counter + 1;
0211                     
0212                     
0213                     
0214                 end
0215                 
0216                 if isSingleRegC && startcounter~=counter % one cell to be mapped
0217                     totCost(:,startcounter:counter-1) = areaChangeFactor * 1./areaOverlapTransCost(:,startcounter:counter-1) + ...
0218                         areaFactor * 1./areaOverlapCost(:,startcounter:counter-1) + centroidCost(:,startcounter:counter-1) + ...
0219                         areaChangeFactor * abs(areaChange(:,startcounter:counter-1));
0220                     
0221                     [~,indx] = min( totCost(:,startcounter:counter-1));
0222                     
0223                     goodOneToOne(ii) = abs(areaChange(startcounter+indx-1)) < 0.15 &&...
0224                         areaOverlapCost(startcounter+indx-1) > 0.7 ...
0225                         && areaOverlapTransCost(startcounter+indx-1) > 0.8;
0226                     
0227                 else
0228                     goodOneToOne (ii) = 1; % no two to two mappings
0229                 end
0230                 
0231                 
0232                 % if the one to one mapping looks good don't bother with
0233                 % pairs - to make faster
0234                 
0235                 if ~goodOneToOne(ii)
0236                     for yy = 1:numel(possibleMapInd)
0237                         for kk = (yy+1) : numel(possibleMapInd)
0238                             
0239                             sis(1) = possibleMapInd(yy);
0240                             sis(2) = possibleMapInd(kk);
0241                             
0242                             isItPair = all(ismember(allF,[sis(1),sis(2)]));
0243                             if any(isItPair)
0244                                 % find their location
0245                                 
0246                                 indexC(:,counter) = cRegs;
0247                                 indexF(:,counter) = sis;
0248                                                                 
0249                                 % combined masks, areas, centroids
0250                                 [maskF,areaF,centroidF] = regProperties (data_f,sis,BB_c_xx,BB_c_yy);
0251                                 overlapMask = maskF(maskC);
0252                                 areaOverlap = sum(overlapMask(:)); % area of overlap between jj and ii
0253                                 areaOverlapCost(counter) = areaOverlap/areaC;
0254                                 
0255                                 if  (areaOverlapCost(counter) == 0)
0256                                     areaOverlapCost(counter) = noOverlap;
0257                                 end
0258                                 
0259                                 displacement = centroidC - centroidF;
0260                                 centroidCost(counter) = sqrt(sum((displacement).^2));
0261                                 
0262                                 if centroidCost == 0
0263                                     outwardMot(counter) = 0;
0264                                 else
0265                                     outwardMot(counter) = (distFromColony*displacement')/centroidCost(counter);
0266                                 end
0267                                 distFromColonyMat (counter) = exp(-distFromColn(ii)/outwardMotFactor);
0268                                 
0269                                 offset = round(displacement);
0270                                 maskOut = imshift(maskF,offset);
0271                                 maskOut = maskOut(maskC);
0272                                 areaOverlapTrns = sum(maskOut(:));
0273                                 areaOverlapTransCost(counter) = areaOverlapTrns/areaC;
0274                                 
0275                                 if  (areaOverlapTransCost(counter) == 0)
0276                                     areaOverlapTransCost(counter) = noOverlap;
0277                                 end
0278                                 
0279                                 areaChange(counter) = (areaF - areaC)/(areaC);
0280                                 counter = counter + 1;
0281                             end
0282                         end
0283                     end
0284                 end
0285             end
0286         end
0287         
0288         
0289         areaChangePenalty = zeros(size(areaChange,1),size(areaChange,2));
0290         
0291         if forward
0292             % area decreases
0293             areaChangePenalty((areaChange) < -0.1) = 100;
0294         else
0295             outwardMot = - outwardMot;
0296             areaChangePenalty((areaChange) > 0.1) = 100;
0297         end
0298         
0299         %  penalty for big area changes
0300         areaChangePenalty(abs(areaChange) > 0.6) = 1000;
0301         areaChangePenalty(abs(areaChange) > 0.3) = 50;
0302         overlapCost = 1 - areaOverlapCost;
0303         totCost = areaChangePenalty +  areaChangeFactor * 1./areaOverlapTransCost + ...
0304             centroidWeight * centroidCost +  areaChangeFactor * abs(areaChange) + ...
0305             distFromColonyMat * areaFactor * 1./areaOverlapCost + outwardMotFactor * outwardMot ;
0306         
0307         %nonNanCost = ~isnan(totCost);
0308         totCost = totCost(1:counter-1);
0309         indexC = indexC(:,1:counter-1);
0310         indexF = indexF(:,1:counter-1);
0311                 
0312         costMat= totCost;        
0313         flagger = ~isnan(costMat);
0314         
0315         while any( flagger(:)) >0
0316             
0317             [~,ind] = min(costMat(:));
0318             assignTemp = indexF(:,ind)';
0319             assignTemp = assignTemp (~isnan(assignTemp));
0320             regionsInC = indexC (:,ind);            
0321             assignments {regionsInC(1)} = assignTemp;
0322             if ~isnan(regionsInC(2))
0323                 assignments {regionsInC(2)} = assignTemp;
0324             end
0325             
0326             %displayMap (data_c,data_f, assignTemp, regionsInC,[],[])
0327             
0328             % find all columns to be set as nans
0329             colToDelF = any(ismember(indexF,assignTemp));
0330             colToDelC = any(ismember(indexC,regionsInC));
0331             costMat (colToDelC) = NaN; % add nans to already assigned
0332             costMat (colToDelF) = NaN; % add nans to already assigned
0333             flagger(colToDelC) = 0;
0334             flagger(colToDelF) = 0;            
0335         end
0336         
0337         % make list of revAssign
0338         revAssign = getRevAssign();
0339         
0340         cArea = [data_c.regs.props.Area];
0341         fArea = [data_f.regs.props.Area];
0342        
0343         % attempt to fix assignment error for cells left without assignment
0344         [assignments,revAssign] =fixProblems(assignments,revAssign, overlapCost, indexF,indexC, cArea, fArea);
0345         [revAssign,assignments] =fixProblems(revAssign,assignments, overlapCost, indexC, indexF,fArea, cArea);        
0346         [assignments,revAssign] = exchangeAssignment (assignments,revAssign, totCost, indexF, indexC);
0347         [revAssign,assignments] = exchangeAssignment (revAssign,assignments, totCost, indexC, indexF);
0348         
0349         dA = changeInArea(assignments, cArea,fArea);
0350         errorR = setError(dA);
0351         
0352         if debug_flag
0353             visualizeLinking(data_c,data_f,assignments);
0354         end
0355         
0356     end
0357 end
0358 
0359 
0360 
0361 
0362 
0363     function [assignments, revAssign] = fixProblems (assignments, revAssign, ...
0364             totCost, indexF, indexC, cArea, fArea)
0365         % fixProblems : used to fix cells not assigned to anything.
0366         
0367         leftInF = find(cellfun('isempty',revAssign));
0368         
0369         for jj = leftInF
0370             
0371             bestAssgnC = findBestSingleAssign (jj, totCost, indexF, indexC) ;
0372             
0373             if ~isnan(bestAssgnC) && bestAssgnC <= numel(assignments)
0374                 FAlready = assignments{bestAssgnC};
0375                 if isempty(FAlready)
0376                     assignments{bestAssgnC} = jj;
0377                 else
0378                     revToAlreadyF = revAssign{FAlready};
0379                     areaC = sum(cArea(revToAlreadyF));
0380                     areaFBefore = sum(fArea(FAlready));
0381                     dABefore = (areaFBefore - areaC)/max(areaFBefore,areaC);
0382                     
0383                     if numel(revToAlreadyF) == 2 && ...
0384                             setError(dABefore)>0
0385                         % two assigned to other f - steal one
0386                         areaFjj = fArea(jj);
0387                         newRevToAlreadyF = revToAlreadyF(revToAlreadyF~=bestAssgnC);
0388                         newAreaC = cArea(newRevToAlreadyF);
0389                         areaC = cArea(bestAssgnC);
0390                         newdAjj = (areaFjj - areaC)/max(areaFjj,areaC);
0391                         newdAalreadyF = (areaFBefore - newAreaC)/max(areaFBefore,areaC);;
0392                         if  ~setError(newdAjj) && ...
0393                                 ~setError(newdAalreadyF)
0394                             assignments{bestAssgnC} = jj;
0395                             revAssign{jj} = bestAssgnC;
0396                             revAssign{FAlready} = newRevToAlreadyF;
0397                         end
0398                     else
0399                         % see if assigning both to bestAssgnC solves the problem
0400                         tempAssgn = [FAlready,jj];
0401                         areaF = areaFBefore + fArea(jj);
0402                         dAtmp = (areaF - areaC)/max(areaF,areaC);
0403                         if  setError(dABefore) > 0 && ...
0404                                 ~setError(dAtmp)
0405                             % disp (['Assign both  ' ,num2str(jj) , '  ', num2str(FAlready), ' to ' , num2str(bestAssgnC)])
0406                             assignments{bestAssgnC} = tempAssgn;
0407                             revAssign{jj} = bestAssgnC;
0408                         end
0409                     end
0410                 end
0411             end
0412             
0413         end
0414         
0415     end
0416 
0417     function bestAssignC = findBestSingleAssign (value_f, totCost, indexF, indexC)
0418         
0419         fAssign =  ((indexF(1,:)== value_f) & (isnan(indexF(2,:))));
0420         totCost(~fAssign) = NaN;
0421         if sum(~isnan(totCost))
0422             [~,indC] = min(totCost);
0423             bestAssignC = indexC (1,indC);
0424         else
0425             bestAssignC = NaN;
0426         end
0427         
0428     end
0429 
0430     function errorR = setError(DA)
0431         if numel(DA) > 0
0432             errorR = zeros(1, numel(DA));
0433             errorR (DA < minDA ) = 2;
0434             errorR (DA > maxDA) = 3;
0435         else
0436             errorR = [];
0437         end
0438     end
0439 
0440     function dA = changeInArea(assignments, cArea,fArea)
0441         % change in area set for data_c
0442         numRegs1 = size(assignments,2);
0443         dA = nan*zeros(1,numRegs1);
0444         for ll = 1 : numRegs1
0445             tmpAssgn =  assignments{ll};
0446             carea_tmp =  (cArea(ll));
0447             farea_tmp = sum(fArea(tmpAssgn));
0448             dA(ll) = (farea_tmp - carea_tmp) / max(carea_tmp,farea_tmp);
0449         end
0450         
0451     end
0452 
0453     function revAssign = getRevAssign()
0454         revAssign = cell( 1, numRegs2);
0455         for ll = 1 : numRegs1
0456             tmpAss =  assignments{ll};
0457             for uu = tmpAss
0458                 revAssign{uu} = [revAssign{uu},ll];
0459             end
0460         end
0461     end
0462 
0463 
0464 
0465     function [assignments,revAssign] = exchangeAssignment (assignments,...
0466             revAssign, totCost, indexF, indexC)
0467         % finds cells without assignment, and finds their best possible assignment.
0468         % exchanges assignments for the one that has taken that assignment if
0469         % it is second best possible choice was also left without assignment.
0470         
0471         leftInC = find(cellfun('isempty',assignments));
0472         leftInF = find(cellfun('isempty',revAssign));
0473         
0474         for kk = 1 : numel(leftInC)
0475             currentCreg = leftInC(kk);
0476             
0477             bestF = findBestSingleAssign (currentCreg, totCost, indexC, indexF);
0478             if ~isempty(bestF)
0479        
0480             for badC = 1 : numel(assignments)
0481                 tempAss = assignments{badC};
0482                 if ~isempty(tempAss) && any(tempAss ==bestF)
0483                     break
0484                 end
0485             end
0486             
0487             flaggerC = (indexC(1,:) == badC) & (isnan(indexC(2,:)));
0488             flaggerF = (indexF(1,:) == bestF) & (isnan(indexF(2,:)));
0489             totCostTemp = totCost;
0490             totCostTemp (~flaggerC) = NaN;
0491             totCostTemp(flaggerF) = NaN;
0492             [~,ind] = min(totCostTemp);
0493             badCSecondF = indexF(:,ind);
0494             
0495             if ~isempty(badCSecondF) && isnan(badCSecondF(2)) &&  any(leftInF == badCSecondF(1))
0496                % disp (['Exchanging assignment for ', num2str(badC), '  and  ', num2str(currentCreg),  ' ', num2str(bestF)]);
0497                 assignments{badC} = badCSecondF(1);
0498                 revAssign{badCSecondF(1)} = badC;
0499                 revAssign{bestF} = currentCreg;
0500                 assignments{currentCreg} = bestF;
0501             end
0502             end
0503         end
0504     end
0505 
0506 
0507 
0508 
0509 end
0510 
0511 
0512 
0513 function [bbx,bby] = getBoxLimits (data_c,regNums)
0514 % getBoxLimits : returns the bounding box for regNums
0515 
0516 regNums = regNums(~isnan(regNums));
0517 comboBoundingBox = [];
0518 ss = size(data_c.phase);
0519 % get total bounding box
0520 for ii = 1: numel(regNums)
0521     reg = regNums(ii);
0522     comboBoundingBox = addBB(comboBoundingBox,data_c.regs.props(reg).BoundingBox);
0523 end
0524 
0525 pad = 20;
0526 [bbx,bby] =  getBBpad(comboBoundingBox,ss,pad);
0527 
0528 end
0529 
0530 
0531 
0532 function [pairsC,neighF] = findNeighborPairs (data_c, numRegs1, regsInC, data_f)
0533 % findNeighborPairs : finds neighboring regions to be considered as pairs
0534 
0535 global str8
0536 counter = 1;
0537 pairsC = NaN * zeros(2,numRegs1 * numRegs1);
0538 neighF = cell( 1, numRegs1);
0539 labels_c = data_c.regs.regs_label;
0540 if ~isempty(data_f)
0541     labels_f = data_f.regs.regs_label;
0542 end
0543 
0544 for jj = regsInC
0545     [bbx,bby] = getBoxLimits (data_c,jj);
0546     labels_c_box = labels_c(bby,bbx);
0547     maskc = (labels_c_box==jj);
0548     tmp_mask = imdilate(maskc, str8);
0549     neigh = labels_c_box(tmp_mask);% get neighbors
0550     
0551     
0552     ind_neigh = unique(neigh)'; % unique neighbors
0553     ind_neigh = ind_neigh(ind_neigh > jj); % not already made pairs
0554     
0555     if  ~isempty(data_f)
0556         labels_f_box = labels_f(bby,bbx);
0557         fneigh = labels_f_box(tmp_mask);
0558         fneigh = fneigh (fneigh~=0);
0559         neighF {jj} = unique(fneigh)';
0560     end
0561     
0562     % make pairs
0563     for uu = ind_neigh
0564         pairsC(:,counter) = [jj;uu];
0565         counter = counter + 1;
0566     end
0567 end
0568 
0569 cleanpairIDs = (nansum(pairsC)~=0);
0570 pairsC = pairsC(:,cleanpairIDs);
0571 
0572 end
0573 
0574 
0575 function [comboMask,comboArea,comboCentroid] = regProperties (data_c,regNums,bbx,bby)
0576 % regProperties :  calculates regNums properties : area, mask and centroid
0577 
0578 comboCentroid = 0;
0579 comboArea = 0;
0580 regNums = regNums(~isnan(regNums));
0581 regs_labels =  data_c.regs.regs_label(bby,bbx);
0582 comboMask = 0 * (regs_labels);
0583 
0584 for ii = 1: numel(regNums)
0585     reg = regNums(ii);
0586     comboCentroid = comboCentroid + data_c.regs.props(reg).Centroid;
0587     comboMask =  comboMask + (regs_labels==reg);
0588     comboArea = comboArea + data_c.regs.props(reg).Area;
0589 end
0590 
0591 comboCentroid = comboCentroid/numel(regNums); % mean centroid
0592 comboMask = (comboMask>0);
0593 
0594 end
0595 
0596 
0597 %         [~,minIndxC] = min(totCost,[],1);
0598 %         % check if it can be added to an assignment?
0599 %         for kk = 1 : numel(leftInF)
0600 %             fToAssign = leftInF(kk);
0601 %             bestAssgn = minIndxC(fToAssign);
0602 %             fAlready = assignments{bestAssgn};
0603 %             tempAssgn = [fAlready,fToAssign];
0604 %             areaF = data_f.regs.props(fAlready).Area + data_f.regs.props(fToAssign).Area;
0605 %             areaC = data_c.regs.props(bestAssgn).Area;
0606 %             dAtmp = (areaF - areaC)/(areaC);
0607 %
0608 %             if  setError(dA(bestAssgn),minDA,maxDA) > 0 && ...
0609 %                  ~setError(dAtmp,minDA,maxDA)
0610 %                 assignments{bestAssgn} = tempAssgn;
0611 %             end
0612 %         end
0613 
0614 
0615 
0616 function visualizeLinking(data_c,data_f,assignments)
0617 figure(1)
0618 clf;
0619 subplot(1,2,1)
0620 imshow(data_c.mask_cell);
0621 subplot(1,2,2)
0622 imshow(data_f.mask_cell);
0623 num_ass = numel(assignments);
0624 randcolor = hsv(256);
0625 markers = {'o','s','d','>','<','^','v','p','h'};
0626 
0627 for c = 1 : num_ass
0628     assF = assignments {c};
0629     if ~isempty(assF)
0630         randomMarker = markers{randi(numel(markers),1)};
0631         randjet = randi(256,1);
0632         color = randcolor(randjet,:);
0633         randjet2 = randi(256,1);
0634         color2 = randcolor(randjet2,:);
0635         figure(1);
0636         subplot(1,2,1)
0637         hold on;
0638         plot(data_c.regs.props(c).Centroid(1),data_c.regs.props(c).Centroid(2),[randomMarker,'k'],'MarkerFaceColor',color,'MarkerEdgeColor',color2,'MarkerSize',8);
0639         subplot(1,2,2)
0640         for i = 1 : numel(assF)
0641             hold on;
0642             plot(data_f.regs.props(assF(i)).Centroid(1),data_f.regs.props(assF(i)).Centroid(2),[randomMarker,'k'],'MarkerFaceColor',color,'MarkerEdgeColor',color2,'MarkerSize',8);
0643         end
0644     end
0645 end
0646 hold on;
0647 subplot(1,2,1)
0648 title('data-c')
0649 subplot(1,2,2)
0650 title('data-f')
0651 
0652 end