errorRezNew

 errorRezNew : links cells from the frame before to the current and attempts to
 resolve segmentation errors if the linking is inconsistent.

 INPUT :
   time : current frame number
   data_c : current time frame data (seg/err) file.
   data_r : reverse time frame data (seg/err) file.
   data_f : forward time frame data (seg/err) file.
   CONST : segmentation parameters.
   cell_count : last cell id used.
   header : last cell id used.
   debug_flag : 1 to display figures for debugging

 OUTPUT :
   data_c : updated current time frame data (seg/err) file.
   data_r : updated reverse time frame data (seg/err) file.
   cell_count : last cell id used.
   resetRegions : if true, regions were modified and this frame needs to
   be relinked.


 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_c, data_r, cell_count,resetRegions] =  errorRezNew (time, ...
0002     data_c, data_r, data_f, CONST, cell_count, header, ignoreError, debug_flag)
0003 % errorRezNew : links cells from the frame before to the current and attempts to
0004 % resolve segmentation errors if the linking is inconsistent.
0005 %
0006 % INPUT :
0007 %   time : current frame number
0008 %   data_c : current time frame data (seg/err) file.
0009 %   data_r : reverse time frame data (seg/err) file.
0010 %   data_f : forward time frame data (seg/err) file.
0011 %   CONST : segmentation parameters.
0012 %   cell_count : last cell id used.
0013 %   header : last cell id used.
0014 %   debug_flag : 1 to display figures for debugging
0015 %
0016 % OUTPUT :
0017 %   data_c : updated current time frame data (seg/err) file.
0018 %   data_r : updated reverse time frame data (seg/err) file.
0019 %   cell_count : last cell id used.
0020 %   resetRegions : if true, regions were modified and this frame needs to
0021 %   be relinked.
0022 %
0023 %
0024 % Copyright (C) 2016 Wiggins Lab
0025 % Written by Stella Stylianidou
0026 % University of Washington, 2016
0027 % This file is part of SuperSegger.
0028 %
0029 % SuperSegger is free software: you can redistribute it and/or modify
0030 % it under the terms of the GNU General Public License as published by
0031 % the Free Software Foundation, either version 3 of the License, or
0032 % (at your option) any later version.
0033 %
0034 % SuperSegger is distributed in the hope that it will be useful,
0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0037 % GNU General Public License for more details.
0038 %
0039 % You should have received a copy of the GNU General Public License
0040 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0041 
0042 global SCORE_LIMIT_MOTHER
0043 global SCORE_LIMIT_DAUGHTER
0044 global REMOVE_STRAY
0045 global header_string
0046 global regToDelete
0047 header_string = header;
0048 verbose = CONST.parallel.verbose;
0049 MIN_LENGTH = 10;
0050 REMOVE_STRAY = CONST.trackOpti.REMOVE_STRAY;
0051 SCORE_LIMIT_DAUGHTER =  CONST.trackOpti.SCORE_LIMIT_DAUGHTER;
0052 SCORE_LIMIT_MOTHER = CONST.trackOpti.SCORE_LIMIT_MOTHER;
0053 DA_MIN = CONST.trackOpti.DA_MIN;
0054 DA_MAX =  CONST.trackOpti.DA_MAX;
0055 regToDelete = [];
0056 resetRegions = false;
0057 
0058 % set all ids to 0
0059 data_c.regs.ID = zeros(1,data_c.regs.num_regs);
0060 modRegions = [];
0061 for regNum =  1 : data_c.regs.num_regs;
0062     
0063     if data_c.regs.ID(regNum) ~= 0
0064         disp ([header, 'ErRes: Frame: ', num2str(time), ' already has an id ',num2str(regNum)]);
0065     elseif ismember (regNum,modRegions)
0066         disp ([header, 'ErRes: Frame: ', num2str(time), ' already modified ',num2str(regNum)]);
0067     else
0068         
0069         rCellsFromC = data_c.regs.map.r{regNum}; % where regNum maps in reverse
0070         
0071         if ~isempty(rCellsFromC)
0072             cCellsTransp = unique([data_c.regs.revmap.r{rCellsFromC}]);
0073             cCellsFromR = unique([data_r.regs.map.f{rCellsFromC}]);
0074         else
0075             cCellsTransp = [];
0076             cCellsFromR = [];
0077         end
0078         
0079         
0080         if ~isempty(cCellsFromR)
0081             rCellsTransp = unique([data_r.regs.revmap.f{cCellsFromR}]);
0082         else
0083             rCellsTransp = [];
0084         end
0085         
0086         
0087         numberMatch = (numel(rCellsFromC) == numel(rCellsTransp)) && ...
0088             (numel(cCellsFromR) == numel(cCellsTransp));
0089         assignmentMatch = numberMatch && all(rCellsFromC == rCellsTransp) && ...
0090             all(cCellsFromR == cCellsTransp);
0091         partialMatch = any(ismember(rCellsFromC,rCellsTransp)) && any(ismember(cCellsFromR, cCellsTransp));
0092         
0093         zeroToOne = numel(rCellsFromC) == 0 ;
0094         oneToOne = numel(rCellsFromC) == 1 &&  numel (cCellsFromR) == 1 ;
0095         oneToTwo = numel(rCellsFromC) == 1 &&  numel (cCellsFromR) == 2 ;
0096         twoToOne = numel(rCellsFromC) == 2 &&  numel (cCellsFromR) == 1 ;
0097         oneToThree = numel(rCellsFromC) == 1 &&  numel (cCellsFromR) == 3 ;
0098         
0099         if numberMatch && assignmentMatch
0100             
0101             if zeroToOne % maps to 0 in the previous frame - stray
0102                 % think whether this is useful :  numel(mapRC) == 0
0103                 if (time ~= 1) && (hasNoFwMapping(data_c, data_f,regNum) && REMOVE_STRAY)
0104                     % deletes the regions not appearing at time = 1 that do not map to anything
0105                     % or if remove_stray flag is set to true.
0106                     data_c.regs.error.label{regNum} = ['Frame: ', num2str(time), ...
0107                         ', reg: ', num2str(regNum), '. is a stray region - Deleted.'];
0108                     if verbose
0109                         disp([header, 'ErRes: ',data_c.regs.error.label{regNum}] );
0110                     end
0111                     regToDelete = [regToDelete;regNum];
0112                     resetRegions = true;
0113                 else % maps to a region in the next frame, or time is 1
0114                     if time~=1
0115                         data_c.regs.error.label{regNum} = ['Frame: ', num2str(time), ...
0116                             ', reg: ', num2str(regNum), '. is a stray region.'];
0117                         if verbose
0118                             disp([header, 'ErRes: ',data_c.regs.error.label{regNum}] );
0119                         end
0120                     end
0121                     [data_c,cell_count] = createNewCell (data_c, regNum, time, cell_count);
0122                 end
0123                 
0124             elseif oneToOne
0125                 % one to one and agreement
0126                 [data_c, data_r] = continueCellLine( data_c, regNum, data_r, rCellsFromC, time, hasError(data_c, regNum));
0127                 
0128                 
0129             elseif oneToTwo
0130                 % one to two : possible splitting event
0131                 
0132                 mother = rCellsFromC;
0133                 sister1 = regNum;
0134                 sister2 = cCellsFromR (cCellsFromR~=regNum);
0135                 haveNoMatch = (isempty(data_c.regs.map.f{sister1}) || isempty(data_c.regs.map.f{sister2}));
0136                 matchToTheSame = ~haveNoMatch && all(ismember(data_c.regs.map.f{sister1}, data_c.regs.map.f{sister2}));
0137                 
0138                 % r: one has no forward mapping, or both map to the same in fw
0139                 if  ~isempty(data_f) && (haveNoMatch || matchToTheSame)
0140                     % wrong division merge cells
0141                     if ~ignoreError
0142                         [data_c,reset_tmp] = merge2Regions (data_c, sister1, sister2, CONST);
0143                         modRegions = [modRegions;sister1;sister2];
0144                         resetRegions = (resetRegions || reset_tmp);
0145                     else
0146                         [data_c,data_r,cell_count,reset_tmp,modids_tmp] = mapBestOfTwo (data_c, cCellsFromR, data_r, rCellsFromC, time, verbose, cell_count,header);
0147                         modRegions = [modRegions;modids_tmp];
0148                     end
0149                     resetRegions = (resetRegions || reset_tmp);
0150                 else
0151                     [data_c, data_r, cell_count] = createDivision (data_c,data_r,mother,sister1,sister2, cell_count, time,header, verbose);
0152                     modRegions = [modRegions;sister1;sister2];
0153                 end
0154                 
0155                 
0156             elseif twoToOne
0157                 % 1 in current maps to two in reverse
0158                 % try to find a segment that should be turned on in current
0159                 
0160                 % The two in reverse map to regNum only : this may be always true by definition
0161                 %twoInRMapToCOnly = numel(data_r.regs.map.f{rCellsFromC(1)}) == 1 && data_r.regs.map.f{rCellsFromC(1)}==regNum && ...
0162                 %   numel(data_r.regs.map.f{rCellsFromC(2)}) == 1 && data_r.regs.map.f{rCellsFromC(2)}==regNum;
0163                 
0164                 success = false;
0165                 if ~ignoreError
0166                     [data_c,success] = missingSeg2to1 (data_c,regNum,data_r,rCellsFromC,CONST);
0167                 end
0168                 
0169                 if success % segment found
0170                     data_c.regs.error.r(regNum) = 0;
0171                     data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0172                         ', reg: ', num2str(regNum),'. Segment added to fix 2 to 1 error'];
0173                     
0174                     if verbose
0175                         disp([header, 'ErRes: ', data_c.regs.error.label{regNum}]);
0176                     end
0177                     if debug_flag
0178                         imshow(cat(3,ag(data_c.regs.regs_label == regNum)+0.5*ag(data_c.phase),...
0179                             ag(data_r.regs.regs_label == rCellsFromC(1)),...
0180                             ag(data_r.regs.regs_label == rCellsFromC(2))));
0181                         keyboard;
0182                     end
0183                     resetRegions = true;
0184                 else
0185                     % maybe copy from frame instead
0186                     % ERROR NOT FIXED : link to the one with the best score
0187                     if debug_flag
0188                     keyboard;
0189                 end
0190                     [data_c,data_r] = mapToBestOfTwo (data_c, regNum, data_r, rCellsFromC, time, verbose,header);
0191                 end
0192             elseif oneToThree
0193                 disp ('merge two?')
0194                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0195                 if debug_flag
0196                     keyboard;
0197                 end
0198             else
0199                 
0200                 disp ('there is another case of correctness?!?')
0201                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0202                 
0203                 if debug_flag
0204                     keyboard;
0205                 end
0206             end
0207         elseif numberMatch
0208             if zeroToOne
0209                 % one to one but disagreement.
0210                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0211                 data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0212                     ', reg: ', num2str(regNum),'. 0->1 Error not fixed - misalign.'];
0213                 
0214                 if debug_flag
0215                     keyboard;
0216                 end
0217             elseif oneToOne
0218                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0219                 data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0220                     ', reg: ', num2str(regNum),'. 1->1 Error not fixed - misalign.'];
0221                 
0222                 if debug_flag
0223                     keyboard;
0224                 end
0225             elseif oneToTwo &&  ~any(cCellsFromR==regNum)
0226                 % regNum maps to mother. but mother maps to two other cells.
0227                 % OTHER POSSIBLE RESOLUTIONS.. :
0228                 % 1 : merging missing, cell divided but piece fell out - check
0229                 % if all three should be mapped
0230                 % 2 : get the best two couples of the three
0231                 
0232                 % create new cell with error
0233                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0234                 [data_c,cell_count] = createNewCell (data_c, regNum, time, cell_count);
0235                 data_c.regs.error.r(regNum) = 1;
0236                 data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0237                     ', reg: ', num2str(regNum),'. Incorrect Mapping 1 to 2 - making a new cell'];
0238                 
0239                 if verbose
0240                     disp([header, 'ErRes: ', data_c.regs.error.label{regNum}]);
0241                 end
0242             else
0243                 disp ('there is another case of matching numbers & un-correctness?!?')
0244                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0245                 
0246                 if debug_flag
0247                     keyboard;
0248                 end
0249             end
0250         else % mistmatch of numbers and assignments
0251             c_matches = cCellsFromR(ismember(cCellsFromR, cCellsTransp));
0252             r_matches = rCellsFromC(ismember(rCellsFromC, rCellsTransp));
0253             
0254             if partialMatch && numel(c_matches) == 1 && numel(r_matches) == 1
0255                 [data_c, data_r] = continueCellLine( data_c, c_matches, data_r, r_matches, time, hasError(data_c, c_matches));
0256                 modRegions = [modRegions;c_matches];
0257                 
0258                 
0259             elseif numel(rCellsFromC) == 1 && numel(cCellsFromR)  == 1 &&  numel (cCellsTransp) == 2
0260                 % c --> r 1, c--->r 1, two c's map to r : inconsistent numbers
0261                 sister1 = regNum;
0262                 sister2 = cCellsTransp (cCellsTransp~=regNum);
0263                 mapRC = cCellsTransp;
0264                 mother = rCellsFromC;
0265                 
0266                 totAreaC = data_c.regs.props(sister1).Area + data_c.regs.props(sister2).Area;
0267                 totAreaR =  data_r.regs.props(mother).Area;
0268                 AreaChange = (totAreaC-totAreaR)/totAreaC;
0269                 divAreaChange = (AreaChange > DA_MIN && AreaChange < DA_MAX);
0270                 haveNoMatch = (isempty(data_c.regs.map.f{sister1}) || isempty(data_c.regs.map.f{sister2}));
0271                 matchToTheSame = ~haveNoMatch && all(ismember(data_c.regs.map.f{sister1}, data_c.regs.map.f{sister2}));
0272                 oneIsSmall = (data_c.regs.info(sister1,1) < MIN_LENGTH) ||  (data_c.regs.info(sister1,1) < MIN_LENGTH);
0273                 
0274                 if divAreaChange && ~ignoreError && ~isempty(data_f) && (haveNoMatch || matchToTheSame || oneIsSmall)
0275                     % r: one has no forward mapping, or both map to the same in fw, or one small
0276                     % wrong division merge cells
0277                     [data_c,mergeReset] = merge2Regions (data_c, sister1, sister2, CONST);
0278                     modRegions = [modRegions;sister1;sister2];
0279                     resetRegions = (resetRegions || mergeReset);
0280                 elseif divAreaChange
0281                     [data_c, data_r, cell_count] = createDivision (data_c,data_r,mother,sister1,sister2, cell_count, time,header, verbose);
0282                     modRegions = [modRegions;sister1;sister2];
0283                 else
0284                     % map to best, remove mapping from second
0285                     [data_c,data_r,cell_count,reset_tmp,modids_tmp] = mapBestOfTwo (data_c, mapRC, data_r, rCellsFromC, time, verbose, cell_count,header);
0286                     resetRegions = (resetRegions || reset_tmp);
0287                     modRegions = [modRegions;modids_tmp];
0288                 end
0289                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0290                 
0291                 if debug_flag
0292                     keyboard;
0293                 end
0294             elseif oneToTwo && any([data_c.regs.map.r{cCellsFromR(cCellsFromR~=regNum)}] ~= rCellsFromC)
0295                 % regNum -> mother, mother maps to two, second does not map to mother.
0296                 [data_c, data_r] = continueCellLine( data_c, regNum, data_r, rCellsFromC, time, hasError(data_c, regNum));
0297                 displayMap (data_c,data_r, rCellsFromC, cCellsTransp,cCellsFromR,rCellsTransp)
0298                 
0299                 if debug_flag
0300                     keyboard;
0301                 end
0302             else
0303                 data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0304                     ', reg: ', num2str(regNum),'. Error not fixed mishmatch and misalign.'];
0305                 
0306                 if debug_flag
0307                     keyboard;
0308                 end
0309                 
0310             end
0311             
0312         end
0313         
0314     end
0315 end
0316 
0317 [data_c] = deleteRegions( data_c,regToDelete);
0318 
0319 end
0320 
0321 
0322 
0323 function intDisplay (data_c,regC,data_f,regF)
0324 % intDisplay : displays linking
0325 % reg : maskF
0326 % green : maskC
0327 % blue : all cell masks  in c
0328 
0329 
0330 maskC = data_c.regs.regs_label*0;
0331 for c = 1 : numel(regC)
0332     if ~isnan(regC(c))
0333         maskC = maskC + (data_c.regs.regs_label == regC(c))>0;
0334     end
0335 end
0336 
0337 if ~isempty (data_f)
0338     maskF = data_f.regs.regs_label*0;
0339     if isempty(regF)
0340         disp('nothing')
0341         imshow (cat(3,0*ag(maskF),ag(maskC),ag(data_c.mask_cell)));
0342     else
0343         for f = 1 : numel(regF)
0344             if ~isnan(regF(f))
0345                 maskF = maskF + (data_f.regs.regs_label == regF(f))>0;
0346             end
0347         end
0348         imshow (cat(3,ag(maskF),ag(maskC),ag(data_c.mask_cell)));
0349     end
0350 end
0351 end
0352 
0353 
0354 
0355 
0356 function [ data_c, data_r, cell_count ] = createDivision (data_c,data_r,mother,sister1,sister2, cell_count, time, header, verbose)
0357 global SCORE_LIMIT_MOTHER
0358 global SCORE_LIMIT_DAUGHTER
0359 
0360 data_c.regs.error.label{sister1} = (['Frame: ', num2str(time),...
0361     ', reg: ', num2str(sister1),' and ', num2str(sister2),' . good cell division from mother reg', num2str(mother),'. [L1,L2,Sc] = [',...
0362     num2str(data_c.regs.L1(sister1),2),', ',num2str(data_c.regs.L2(sister1),2),...
0363     ', ',num2str(data_c.regs.scoreRaw(sister1),2),'].']);
0364 if verbose
0365     disp([header, 'ErRes: ', data_c.regs.error.label{sister1}] );
0366 end
0367 data_r.regs.error.r(mother) = 0;
0368 data_c.regs.error.r(sister1) = 0;
0369 data_c.regs.error.r(sister2) = 0;
0370 [data_c, data_r, cell_count] = markDivisionEvent( ...
0371     data_c, sister1, data_r, mother, time, 0, sister2, cell_count);
0372 end
0373 
0374 
0375 function result = hasNoFwMapping (data_c,data_f,regNum)
0376 result = ~isempty(data_f) && isempty(data_c.regs.map.f{regNum});
0377 end
0378 
0379 
0380 function [data_c,data_r] = mapToBestOfTwo (data_c, regNum, data_r, mapCR, time, verbose,header)
0381 % maps to best from two forward
0382 
0383 cost1 = data_c.regs.cost.r(regNum,mapCR(1));
0384 cost2 = data_c.regs.cost.r(regNum,mapCR(2));
0385 
0386 if cost1<cost2 || isnan(cost2)
0387     data_c.regs.error.r(regNum) = 2;
0388     errorStat = 1;
0389     data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0390         ', reg: ', num2str(regNum),'. Map to minimum cost'];
0391     
0392     if verbose
0393         disp([header, 'ErRes: ', data_c.regs.error.label{regNum}]);
0394     end
0395     [data_c, data_r] = continueCellLine(data_c, regNum, data_r, mapCR(1), time, errorStat);
0396     
0397 else
0398     errorStat = 1;
0399     data_c.regs.error.r(regNum) = 2;
0400     data_c.regs.error.label{regNum} = ['Frame: ', num2str(time),...
0401         ', reg: ', num2str(regNum),'. Map to minimum cost'];
0402     
0403     if verbose
0404         disp([header, 'ErRes: ', data_c.regs.error.label{regNum}]);
0405     end
0406     [data_c, data_r] = continueCellLine(data_c, regNum, data_r, mapCR(2), time, errorStat);
0407     
0408 end
0409 end
0410 
0411 
0412 function [data_c,data_r,cell_count,resetRegions,idsOfModRegions] = mapBestOfTwo ...
0413     (data_c, mapRC, data_r, mapCR, time, verbose, cell_count,header)
0414 % maps to best from two forward
0415 global regToDelete
0416 global REMOVE_STRAY
0417 resetRegions = 0;
0418 [~,minInd] = min (data_c.regs.dA.r(mapRC));
0419 keeper = mapRC(minInd);
0420 remove = mapRC(mapRC~=keeper);
0421 [data_c, data_r] = continueCellLine( data_c, keeper, data_r, mapCR, time, hasError(data_c, keeper));
0422 data_c.regs.revmap.r{mapCR} = keeper;
0423 
0424 
0425 data_c.regs.error.r(remove) = 1;
0426 idsOfModRegions = [remove;keeper];
0427 if REMOVE_STRAY && hasNoFwMapping(data_c, data_f,remove)
0428     data_c.regs.error.label{remove} = (['Frame: ', num2str(time),...
0429         ', reg: ', num2str(remove),' was not the best match for ', num2str(mapCR),' and was deleted.' num2str(keeper) , ' was.']);
0430     if verbose
0431         disp([header, 'ErRes: ', data_c.regs.error.label{remove}] );
0432     end
0433     regToDelete = [regToDelete;remove];
0434     resetRegions = true;
0435 else
0436     data_c.regs.error.label{remove} = (['Frame: ', num2str(time),...
0437         ', reg: ', num2str(remove),' was not the best match for ', num2str(mapCR),' .']);
0438     if verbose
0439         disp([header, 'ErRes: ', data_c.regs.error.label{remove}] );
0440     end
0441 end
0442 
0443 end
0444 
0445 function isError = hasError(data_c, regNum)
0446 isError = 0; % FIX ME
0447 %isError = data_c.regs.error.r(regNum) ~= 0 || data_c.regs.error.f(regNum) ~= 0;
0448 end