diff --git a/Internal/intFixPos.m b/Internal/intFixPos.m new file mode 100644 index 0000000..360c45e --- /dev/null +++ b/Internal/intFixPos.m @@ -0,0 +1,17 @@ +function x = intFixPos( x, ss ) + + +if x(1) < 1 + x(1) = 1; +elseif x(1) > ss(2) + x(1) = ss(2); +end + +if x(2) < 1 + x(2) = 1; +elseif x(2) > ss(1) + x(2) = ss(1); +end + + +end diff --git a/cell/trackOptiClist.m b/cell/trackOptiClist.m index 878abb0..8fa73cf 100644 --- a/cell/trackOptiClist.m +++ b/cell/trackOptiClist.m @@ -89,6 +89,9 @@ index_dlminOld = strcmp(tmpFields,'dlmin'); index_ehist = find(strcmp(tmpFields,'error_frame')); + + prog = []; + prog0 = []; % loop through all the images (*err.mat files) for i = 1:num_im data_c = loaderInternal([dirname,contents(i).name]); @@ -182,6 +185,57 @@ daughter1_id = drill(data_c.regs.daughterID,'(1)'); daughter2_id = drill(data_c.regs.daughterID,'(2)'); + mother_id = data_c.regs.motherID; + + % get lineage stuff done here + + if i == 1 + prog = ID; + gen0 = 0*ID; + gen = 0*ID; + + else + % copy info from previous IDs + prog = nan( size( ID ) ); + gen0 = nan( size( ID ) ); + gen = nan( size( ID ) ); + + prog(~birthID_flag) = prog_(ID(~birthID_flag)); + gen0(~birthID_flag) = gen0_(ID(~birthID_flag)); + gen(~birthID_flag) = gen_(ID(~birthID_flag)); + + % now take care of new cells + mother_list = mother_id(birthID_flag); + + prog_tmp = nan( size( mother_list )); + gen0_tmp = nan( size( mother_list )); + gen_tmp = nan( size( mother_list )); + + prog_tmp(mother_list>0) = prog_(mother_list(mother_list>0)); + gen0_tmp(mother_list>0) = gen0_(mother_list(mother_list>0))+1; + gen_tmp(mother_list>0) = gen_(mother_list(mother_list>0))+1; + + ID_ml = ID(birthID_flag); + prog_tmp(mother_list==0) = ID_ml(mother_list==0); + gen0_tmp(mother_list==0) = nan; + gen_tmp(mother_list==0) = 0; + + prog(birthID_flag) = prog_tmp; + gen0(birthID_flag) = gen0_tmp; + gen(birthID_flag) = gen_tmp; + end + + IDmax = max(ID); + + gen0_ = nan( [1,IDmax] ); + gen0_(ID) = gen0; + + gen_ = nan( [1,IDmax] ); + gen_(ID) = gen; + + prog_ = nan( [1,IDmax] ); + prog_(ID) = prog; + % done lineage stuff. length1 = drill(data_c.CellA,'.length(1)'); length2 = drill(data_c.CellA,'.length(2)'); @@ -421,6 +475,9 @@ {'Neck width'},{'neckWidth'},0,0; {'Maximum width'},{'maxWidth'},0,0; {'Cell dist to edge'},{'cell_dist'},1,0; + {'Generation'},{'gen'},0,0; + {'Generation0'},{'gen0'},0,0; + {'Progenitor'},{'prog'},0,0; ]; diff --git a/segmentation/defineGoodSegs.m b/segmentation/defineGoodSegs.m index 2380a17..8183da4 100644 --- a/segmentation/defineGoodSegs.m +++ b/segmentation/defineGoodSegs.m @@ -1,5 +1,4 @@ -function [data] = defineGoodSegs(data, ws, phaseNorm, C2phaseThresh, ... - mask_bg, A, CONST, calcScores) +function [data] = defineGoodSegs(data, ws, CONST, calcScores) % defineGoodSegs sets the segments to good, bad and 3n set by the watershed algorithm % "Good" segments (segs_good) are the ones that lie along a real cellular % boundary, "bad" segments, lie along spurious boundaries @@ -11,9 +10,6 @@ % INPUT : % data : data segmentation frame file % ws : watersehd image -% phaseNorm : normalized phase image -% C2phaseThresh : c2 (curvature) calculated image -% mask_bg : background mask % A : scoring coefficients % CONST : segmentation parameters % caclScores : boolean to calculate scores @@ -39,6 +35,13 @@ % along with SuperSegger. If not, see . +% phaseNorm : normalized phase image +% C2phaseThresh : c2 (curvature) calculated image +% mask_bg : background mask +mask_bg = data.mask_bg; +phaseNorm = data.phaseNorm; +C2phaseThresh = data.C2phaseThresh; +A = CONST.superSeggerOpti.A; sim = size( phaseNorm ); diff --git a/segmentation/intRemovePillars.m b/segmentation/intRemovePillars.m new file mode 100644 index 0000000..727fa67 --- /dev/null +++ b/segmentation/intRemovePillars.m @@ -0,0 +1,63 @@ +function mask = intRemovePillars(phase, CONST) + + + +if exist( 'CONST' ) && ~isempty( CONST ) && ... + isfield( CONST.superSeggerOpti, 'remove_pillars' ) + radius = CONST.superSeggerOpti.remove_pillars.radius; + cut = CONST.superSeggerOpti.remove_pillars.cut; + Area_Cut = CONST.superSeggerOpti.remove_pillars.Area_Cut; + debug = CONST.superSeggerOpti.remove_pillars.debug; + +else + radius = 2; + cut = 0.05; + Area_Cut = 700; + debug = false; +end + + +[~,~,~,D] = curveFilter( double(phase), radius ); + +D = D/max(D(:)); + + +D(D (mult_max*mean_phase)) = mult_max*mean_phase; phaseNorm(phaseNorm < (mult_min*mean_phase)) = mult_min*mean_phase; +% Filter to remove internal structure from phase image ofcells here +% (and little pieces of debris from the background +if isfield( CONST.superSeggerOpti, 'remove_int_struct' ) && ... + CONST.superSeggerOpti.remove_int_struct.flag + tmp_im = imclose( phaseNorm, strel('disk',... + CONST.superSeggerOpti.remove_int_struct.radius )); + W = CONST.superSeggerOpti.remove_int_struct.weight; + phaseNorm = (1-W)*phaseNorm + W*tmp_im; +end + + % if the size of the matrix is even, we get a half pixel shift in the % position of the mask which turns out to be a problem later. @@ -212,6 +225,14 @@ adapt_flag=1; end + +% Remove CellASIC pillars here +if isfield( CONST.superSeggerOpti, 'remove_pillars' ) && ... + CONST.superSeggerOpti.remove_pillars.flag + mask_to_remove = ~intRemovePillars(phaseOrig,CONST); + mask_bg(mask_to_remove) = false; +end + %% Split up the micro colonies into watershed regions to assemble cells % if data exists, reconstruct ws if dataFlag @@ -287,7 +308,12 @@ %% Remake the data structure % Determine the "good" and "bad" segments -data = defineGoodSegs(data,ws,phaseNormUnfilt,C2phaseThresh,mask_bg, A,CONST, ~dataFlag ); +data.mask_bg = mask_bg; +data.phaseNorm = phaseNormUnfilt; +data.C2phaseThresh = C2phaseThresh; + +data = defineGoodSegs(data, ws, CONST, ~dataFlag ); + data.mask_colonies = mask_colonies; % copy the existing score into the data structure @@ -315,7 +341,7 @@ if disp_flag - figure(1) + figure( 'name', 'SuperSegger frame segmentation' ) clf; showSegDataPhase(data); drawnow; diff --git a/trainingConstants/updateTrainingImage.m b/trainingConstants/updateTrainingImage.m index 46e9ae4..5aca34d 100644 --- a/trainingConstants/updateTrainingImage.m +++ b/trainingConstants/updateTrainingImage.m @@ -33,82 +33,154 @@ ss = size(data.phase); x = round(x); +%% +x = intFixPos(x,ss); -if ~isempty(x) - % creates an image of 51 x 51 of gaussian like point - tmp = zeros([51,51]); - tmp(26,26) = 1; - tmp = 8000-double(bwdist(tmp)); - - rmin = max([1,x(2)-25]); - rmax = min([ss(1),x(2)+25]); - - cmin = max([1,x(1)-25]); - cmax = min([ss(2),x(1)+25]); - - rrind = rmin:rmax; - ccind = cmin:cmax; - - pointSize = [numel(rrind),numel(ccind)]; - - - if im_flag == 1 - segs = data.segs.segs_good(rrind,ccind) + ... - data.segs.segs_bad(rrind,ccind); - segs = segs>0; - tmp = tmp(26-x(2)+rrind,26-x(1)+ccind).*segs ; +if ~isempty(x) + + if FLAGS.whichButton == FLAGS.EdgeToggleRadioButtonFlag; + %% in toggle mode + + % creates an image of 51 x 51 of gaussian like point + tmp = zeros([51,51]); + tmp(26,26) = 1; + tmp = 8000-double(bwdist(tmp)); + + rmin = max([1,x(2)-25]); + rmax = min([ss(1),x(2)+25]); + + cmin = max([1,x(1)-25]); + cmax = min([ss(2),x(1)+25]); + + rrind = rmin:rmax; + ccind = cmin:cmax; + + pointSize = [numel(rrind),numel(ccind)]; + + + if im_flag == 1 + + segs = data.segs.segs_good(rrind,ccind) + ... + data.segs.segs_bad(rrind,ccind); + segs = segs>0; + tmp = tmp(26-x(2)+rrind,26-x(1)+ccind).*segs ; + + [~,ind] = max( tmp(:) ); + + % indices in point image for max / closest segment + [sub1, sub2] = ind2sub( pointSize, ind ); + + % closest segments id + ii = data.segs.segs_label(rmin+sub1-1,cmin+sub2-1); + + if ii ~=0 + + % xx and yy are the segments coordinates + [xx,yy] = getBB( data.segs.props(ii).BoundingBox ); + + if data.segs.score(ii) % score is 1 + data.segs.score(ii) = 0; % set to 0 + data.segs.segs_good(yy,xx) ... + = double(~~(data.segs.segs_good(yy,xx)... + - double(data.segs.segs_label(yy,xx)==ii))); + data.segs.segs_bad(yy,xx) = ... + double(~~(data.segs.segs_bad(yy,xx)... + +double(data.segs.segs_label(yy,xx)==ii))); + else + data.segs.score(ii) = 1; + data.segs.segs_good(yy,xx) = ... + double(~~(data.segs.segs_good(yy,xx)+... + double(data.segs.segs_label(yy,xx)==ii))); + data.segs.segs_bad(yy,xx) = ... + double(~~(data.segs.segs_bad(yy,xx)-... + double(data.segs.segs_label(yy,xx)==ii))); + end + + % updates cell mask + data.mask_cell = double((data.mask_bg - data.segs.segs_good - data.segs.segs_3n)>0); + data.regs.regs_label = bwlabel(data.mask_cell); + touch_list = [touch_list, ii]; + end + elseif im_flag == 2 + tmp = tmp(26-x(2)+rrind,26-x(1)+ccind).*data.mask_cell(rrind,ccind); + try + [~,ind] = max( tmp(:) ); + catch ME + printError(ME); + end + + [sub1, sub2] = ind2sub( pointSize, ind ); + ii = data.regs.regs_label(sub1-1+rmin,sub2-1+cmin); + plot( sub2-1+cmin, sub1-1+rmin, 'g.' ); + + if ii + data.regs.score(ii) = ~data.regs.score(ii); + end + end + + + elseif FLAGS.whichButton == FLAGS.RemoveCellRadioButtonFlag + %% Remove cell option + + ii = data.regs.regs_label(x(2),x(1)); + + if ii + data.mask_bg(data.regs.regs_label==ii) = 0; + touch_list = [touch_list, ii]; + + end + intDoUpdate + + + elseif FLAGS.whichButton == FLAGS.BackgroundRadioButtonFlag + %% Background option + + data.mask_bg(x(2),x(1)) = false; + intDoUpdate + + elseif FLAGS.whichButton == FLAGS.CellMaskRadioButtonFlag + %% Cell Mask option + + data.mask_bg(x(2),x(1)) = true; + intDoUpdate + end + + +end + +% updates cell mask + function intDoUpdate + ws = data.segs.segs_good + data.segs.segs_3n + data.segs.segs_bad; + data_p = defineGoodSegs(data, ws, FLAGS.CONST, true ); + + data = intMapEdges( data_p, data ); + + data.mask_cell = double((data.mask_bg - data.segs.segs_good - data.segs.segs_3n)>0); + data.regs.regs_label = bwlabel(data.mask_cell); + + + end +end - [~,ind] = max( tmp(:) ); - % indices in point image for max / closest segment - [sub1, sub2] = ind2sub( pointSize, ind ); - % closest segments id - ii = data.segs.segs_label(rmin+sub1-1,cmin+sub2-1); - if ii ~=0 +function data = intMapEdges( data_p, data ) - % xx and yy are the segments coordinates - [xx,yy] = getBB( data.segs.props(ii).BoundingBox ); - if data.segs.score(ii) % score is 1 - data.segs.score(ii) = 0; % set to 0 - data.segs.segs_good(yy,xx) ... - = double(~~(data.segs.segs_good(yy,xx)... - - double(data.segs.segs_label(yy,xx)==ii))); - data.segs.segs_bad(yy,xx) = ... - double(~~(data.segs.segs_bad(yy,xx)... - +double(data.segs.segs_label(yy,xx)==ii))); - else - data.segs.score(ii) = 1; - data.segs.segs_good(yy,xx) = ... - double(~~(data.segs.segs_good(yy,xx)+... - double(data.segs.segs_label(yy,xx)==ii))); - data.segs.segs_bad(yy,xx) = ... - double(~~(data.segs.segs_bad(yy,xx)-... - double(data.segs.segs_label(yy,xx)==ii))); - end +ns = numel( data_p.segs.score ); - % updates cell mask - data.mask_cell = double((data.mask_bg - data.segs.segs_good - data.segs.segs_3n)>0); - data.regs.regs_label = bwlabel(data.mask_cell); - touch_list = [touch_list, ii]; - end - elseif im_flag == 2 - tmp = tmp(26-x(2)+rrind,26-x(1)+ccind).*data.mask_cell(rrind,ccind); - try - [~,ind] = max( tmp(:) ); - catch ME - printError(ME); - end +for ii = 1:ns + + ind = unique(data.segs.segs_label( data_p.segs.segs_label==ii )); + ind = ind(logical(ind)); + + if ~isempty( ind ) + data_p.segs.score(ii) = data.segs.score(ind(1)); + end +end - [sub1, sub2] = ind2sub( pointSize, ind ); - ii = data.regs.regs_label(sub1-1+rmin,sub2-1+cmin); - plot( sub2-1+cmin, sub1-1+rmin, 'g.' ); +data = data_p; - if ii - data.regs.score(ii) = ~data.regs.score(ii); - end - end end \ No newline at end of file diff --git a/viz/editSegmentsGui.fig b/viz/editSegmentsGui.fig index 2971dba..bfa79f5 100644 Binary files a/viz/editSegmentsGui.fig and b/viz/editSegmentsGui.fig differ diff --git a/viz/editSegmentsGui.m b/viz/editSegmentsGui.m index ccf2bf0..650ba90 100644 --- a/viz/editSegmentsGui.m +++ b/viz/editSegmentsGui.m @@ -26,6 +26,12 @@ 'gui_OutputFcn', @editSegmentsGui_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); + + + + + + if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end @@ -39,13 +45,39 @@ % Global functions -function clickOnImage(hObject, eventdata, handles) -global settings -FLAGS.im_flag = settings.handles.im_flag; -currentData = load([settings.handles.dirname, settings.handles.contents(str2double(settings.handles.frame_no.String)).name]); -[data, list] = updateTrainingImage(currentData, FLAGS, eventdata.IntersectionPoint(1:2)); -save([settings.handles.dirname, settings.handles.contents(str2double(settings.handles.frame_no.String)).name], '-STRUCT', 'data'); -updateUI(settings.hObject, settings.handles); +function clickOnImage(hObject, eventdata, handles, hObject2) +%global settings +disp( 'clickOnImage' ); +hObject +handles + +FLAGS.im_flag = handles.im_flag; +FLAGS.EdgeToggleRadioButtonFlag = 1; +FLAGS.RemoveCellRadioButtonFlag = 2; +FLAGS.BackgroundRadioButtonFlag = 3; +FLAGS.CellMaskRadioButtonFlag = 4; +FLAGS.CONST = handles.CONST; + + +%currentData = load([settings.handles.dirname, settings.handles.contents(str2double(settings.handles.frame_no.String)).name]); + +if handles.EdgeToggleRadioButton.Value + FLAGS.whichButton = FLAGS.EdgeToggleRadioButtonFlag; +elseif handles.RemoveCellRadioButton.Value + FLAGS.whichButton = FLAGS.RemoveCellRadioButtonFlag; +elseif handles.BackgroundRadioButton.Value + FLAGS.whichButton = FLAGS.BackgroundRadioButtonFlag; +elseif handles.CellMaskRadioButton.Value + FLAGS.whichButton = FLAGS.CellMaskRadioButtonFlag; +else + disp( 'Error in radio buttons for editSegments' ); +end + + + +[handles.data, list] = updateTrainingImage( handles.data, FLAGS, eventdata.IntersectionPoint(1:2)); + +updateUI(hObject2, handles); function data = loaderInternal(filename) data = load(filename); @@ -62,17 +94,32 @@ function editSegmentsGui_OpeningFcn(hObject, eventdata, handles, varargin) handles.num_im = length(handles.contents); handles.im_flag = 1; axis tight; + +handles.EdgeToggleRadioButton.Value = 1; +handles.RemoveCellRadioButton.Value = 0; +handles.BackgroundRadioButton.Value = 0; +handles.CellMaskRadioButton.Value = 0; + +handles.data = loaderInternal([handles.dirname, ... + handles.contents(str2double(handles.frame_no.String)).name]); + updateUI(hObject, handles); function updateUI(hObject, handles) -global settings +disp('updateUI'); +hObject +handles + +%global settings delete(get(handles.axes1, 'Children')); -data = loaderInternal([handles.dirname, handles.contents(str2double(handles.frame_no.String)).name]); -data.mask_cell = double((data.mask_bg - data.segs.segs_good - data.segs.segs_3n) > 0); + +data = handles.data; + +handles.data.mask_cell = double((data.mask_bg - data.segs.segs_good - data.segs.segs_3n) > 0); showSegData(data, handles.im_flag, handles.axes1); -settings.handles = handles; -settings.hObject = hObject; -set(handles.axes1.Children, 'ButtonDownFcn', @clickOnImage); +%settings.handles = handles; +%settings.hObject = hObject; +set(handles.axes1.Children, 'ButtonDownFcn', {@clickOnImage, handles, hObject}); guidata(hObject, handles); function figure1_CloseRequestFcn(hObject, eventdata, handles) @@ -91,7 +138,12 @@ function frame_no_Callback(hObject, eventdata, handles) else handles.frame_no.String = num2str(c); end + +handles.data = loaderInternal([handles.dirname, ... + handles.contents(str2double(handles.frame_no.String)).name]); + updateUI(hObject, handles) +guidata(hObject, handles); function frame_no_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) @@ -145,4 +197,96 @@ function relink_Callback(hObject, eventdata, handles) header = 'Relink: '; trackOpti(handles.dirname_xy,skip, handles.CONST, header, startEnd); -end \ No newline at end of file +end + + +% --- Executes on button press in EdgeToggleRadioButton. +function EdgeToggleRadioButton_Callback(hObject, eventdata, handles) +% hObject handle to EdgeToggleRadioButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of EdgeToggleRadioButton + +handles.RemoveCellRadioButton.Value = 0; +handles.BackgroundRadioButton.Value = 0; +handles.CellMaskRadioButton.Value = 0; + +disp('EdgeToggleRadioButton_Callback'); +hObject +handles + +% --- Executes on button press in RemoveCellRadioButton. +function RemoveCellRadioButton_Callback(hObject, eventdata, handles) +% hObject handle to RemoveCellRadioButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of RemoveCellRadioButton + +handles.EdgeToggleRadioButton.Value = 0; +handles.BackgroundRadioButton.Value = 0; +handles.CellMaskRadioButton.Value = 0; + +disp('RemoveCellRadioButton_Callback'); +hObject +handles + +% --- Executes on button press in BackgroundRadioButton. +function BackgroundRadioButton_Callback(hObject, eventdata, handles) +% hObject handle to BackgroundRadioButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of BackgroundRadioButton + +handles.EdgeToggleRadioButton.Value = 0; +handles.RemoveCellRadioButton.Value = 0; +handles.CellMaskRadioButton.Value = 0; + +disp('BackgroundRadioButton_Callback'); +hObject +handles + +% --- Executes on button press in CellMaskRadioButton. +function CellMaskRadioButton_Callback(hObject, eventdata, handles) +% hObject handle to CellMaskRadioButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of CellMaskRadioButton + + +handles.RemoveCellRadioButton.Value = 0; +handles.BackgroundRadioButton.Value = 0; +handles.EdgeToggleRadioButton.Value = 0; + +disp('CellMaskRadioButton_Callback'); +hObject +handles + +% --- Executes on button press in SavePushButton. +function SavePushButton_Callback(hObject, eventdata, handles) +% hObject handle to SavePushButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes on button press in RevertPushButton. + + +data = handles.data; + +save([handles.dirname, handles.contents(... + str2double(handles.frame_no.String)).name], '-STRUCT', 'data'); + + +function RevertPushButton_Callback(hObject, eventdata, handles) +% hObject handle to RevertPushButton (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +handles.data = loaderInternal([handles.dirname, handles.contents(... + str2double(handles.frame_no.String)).name]); + +updateUI(hObject, handles); +guidata(hObject, handles); diff --git a/viz/showSegData.m b/viz/showSegData.m index baa353d..542374e 100644 --- a/viz/showSegData.m +++ b/viz/showSegData.m @@ -29,6 +29,8 @@ function showSegData( data, im_flag, gui_fig ) segs_good = data.segs.segs_good; segs_bad = data.segs.segs_bad; mask_bg = data.mask_bg; +ss = size(mask_bg); + if ~exist('im_flag','var') im_flag = 1; @@ -39,16 +41,30 @@ function showSegData( data, im_flag, gui_fig ) if ~exist('gui_fig','var') || isempty(gui_fig) figure(1); + tmp = [1 1 1 1]; else axes(gui_fig); % new, ie. used in the GUI version + tmp = axis; end if im_flag == 1 % displays good, 3n and bad segments phaseBackag = (ag((~data.mask_cell))); - imshow( cat(3, 0.2*(phaseBackag) + 0.3*ag(segs_3n) + ag(segs_good), ... - 0.2*(phaseBackag) + 0.3*ag(segs_3n) , ... - 0.2*(phaseBackag) + 0.3*ag(segs_3n) + ag(segs_bad) ), 'InitialMagnification', 'fit'); + imshow( cat(3, 0.2*(phaseBackag) + 0.3*ag(segs_3n) + ag(segs_good), ... + 0.2*(phaseBackag) + 0.3*ag(segs_3n) , ... + 0.2*(phaseBackag) + 0.3*ag(segs_3n) + ag(segs_bad) ), ... + 'InitialMagnification', 'fit',... + 'Parent', gui_fig); + +axes( gui_fig ); + +if all( tmp == [0 1 0 1] ); +axis( [1,ss(2),1,ss(1)] ); +else + axis(tmp ) +end + + % 'hi' elseif im_flag == 2 % displays cell mask cc = bwconncomp(cell_mask, 4); labeled = labelmatrix(cc);