Identifying regions in matrix rows

I have an M x N matrix that contains the indices of columns where an external "validity criterion" is fulfilled, and -1 otherwise. Example:

matrix = zeros (5,8) -1;
matrix(1,1:3) = 1:3;
matrix(3,2:3)  = 2:3;
matrix(3,5:7) = 5:7;
matrix(4,2:5) = 2:5;
matrix

matrix = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

Let's call the input above the "regions", i.e., here we have in total 4 regions (always considered row-wise). My objective is to identify, for each row, the largest region, and within that region the central value (rounded up). Currently, I do this as follows:

dummy = zeros(5,1) -1;
jumps = abs(diff([dummy, matrix, dummy], [], 2))

jumps = 5×9

2 1 1 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 4 6 1 1 8 0 0 3 1 1 1 6 0 0 0 0 0 0 0 0 0 0 0 0

[borders_row, borders_col] = find(jumps>1);

That way, I identify the borders of the regions. Afterwards, I go through each row individually to identify the largest region, as I mentioned, which I do as follows:

first_idx = zeros(1,5);
second_idx = false(1,5);
for i = 1:5
    borders_cur = borders_col(borders_row==i).';
    ranges = diff(borders_cur);
    ranges = ranges(1:2:end);
    [width, path_idx] = max(ranges);
    idxx = borders_cur(path_idx)-1 + ceil(width/2);
    if ~isempty(idxx)
        first_idx(i) = idxx;
        second_idx(i) = logical(idxx);
    end
end

For each row, I consider the current borders and first calculate the distance between them, so that I can judge which region is larger if there is more than one. I skip every second value because those indicate the "invalid spaces" in between the regions. I then pick the bigger region and calculate its central index, which I then record in two arrays (first_idx and second_idx).

I want to get rid of the loop as it's computationally expensive, and generally optimise the code as it's part of a bigger loop, and this is one of the slowest parts of the code currently. Any suggestions?

[Matt J]

I want to get rid of the loop as it's computationally expensive

Is it? It only runs from 1 to 5.

[Dominik Rhiem]

Sorry, I should have been more clear. The above is an example that is supposed to showcase the problem. In reality, the loop can easily run from i = 1:20,000. This is also a nested loop, so it will usually be repeated a few thousand times more, at least.

[Bruno Luong]

IMO there is a bug in these three lines

ranges = ranges(1:2:end);
[width, path_idx] = max(ranges);
idxx = borders_cur(path_idx)-1 + ceil(width/2)

Your range is step 2, your path_idx is relative to then subarray of step2 so you should index with borders_cur(2*path_idx-1) in the last line.

Or alternatively you have to make borders_cur step 2 too.

ranges = ranges(1:2:end);
borders_cur = borders_cur(1:2:end); % missing
[width, path_idx] = max(ranges);
idxx = borders_cur(path_idx)-1 + ceil(width/2)

[Dominik Rhiem]

@Bruno Luong Good catch! I have decided to use borders_cur(2*path_idx-1) since it is more clear, in my opinion.

Answer 1

[Bruno Luong]

matrix = zeros (5,8) -1;
matrix(1,1:3) = 1:3;
matrix(3,2:3)  = 2:3;
matrix(3,5:7) = 5:7;
matrix(4,2:5) = 2:5;
matrix

matrix = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

dummy = zeros(5,1) -1;
jumps = abs(diff([dummy, matrix, dummy], [], 2));
[borders_row, borders_col] = find(jumps>1);
first_idx = zeros(1,5);
second_idx = false(1,5);
for i = 1:5
    borders_cur = borders_col(borders_row==i).';
    ranges = diff(borders_cur);
    ranges = ranges(1:2:end);
    [width, path_idx] = max(ranges);
    idxx = borders_cur(2*path_idx-1)-1 + ceil(width/2); % modified by BLU
    if ~isempty(idxx)
        first_idx(i) = idxx;
        second_idx(i) = logical(idxx);
    end
end
first_idx

first_idx = 1×5

2 0 6 3 0

% my method
A = (matrix>0).';
n = size(A,2);
f = false(1, n);
[i,j] = find(diff([f; A; f]));
c = i(1:2:end);
w = i(2:2:end)-c;
r = j(1:2:end);
[rw,is] = sortrows([r,w],[1 -2]);
keep = [true; diff(rw(:,1))>0];
rw = rw(keep,:);
c = c(is(keep));
r = rw(:,1);
w = rw(:,2);
mididx = c-1+ceil(w/2);
first_idx = accumarray(r, mididx, [n,1])'

first_idx = 1×5

2 0 6 3 0

second_idx = logical(first_idx);

[Bruno Luong]

Test on big case

matrix = repmat(1:1000, 10000, 1);
matrix(randi(end, 1000000, 1)) = -1;

tic
n = size(matrix,1);
dummy = zeros(n,1) -1;
jumps = abs(diff([dummy, matrix, dummy], [], 2));
[borders_row, borders_col] = find(jumps>1);
first_idx = zeros(1,n);
second_idx = false(1,n);
for i = 1:n
    borders_cur = borders_col(borders_row==i).';
    ranges = diff(borders_cur);
    ranges = ranges(1:2:end);
    borders_cur = borders_cur(1:2:end); % missing
    [width, path_idx] = max(ranges);
    idxx = borders_cur(path_idx)-1 + ceil(width/2);
    if ~isempty(idxx)
        first_idx(i) = idxx;
        second_idx(i) = logical(idxx);
    end
end
toc

Elapsed time is 32.274977 seconds.

% my method
tic
A = (matrix>0).';
n = size(A,2);
f = false(1, n);
[i,j] = find(diff([f; A; f]));
c = i(1:2:end);
w = i(2:2:end)-c;
r = j(1:2:end);
[rw,is] = sortrows([r,w],[1 -2]);
keep = [true; diff(rw(:,1))>0];
rw = rw(keep,:);
c = c(is(keep));
r = rw(:,1);
w = rw(:,2);
first_idx2 = accumarray(r, c-1+ceil(w/2), [n,1])';
second_idx2 = logical(first_idx2);
toc

Elapsed time is 0.251036 seconds.

isequal(first_idx, first_idx2)

ans = logical
   1

[Dominik Rhiem]

I have been testing your suggestion, and I really like it. It's quick and produces the correct results. Thank you!

[Bruno Luong]

Good ! Glad to help

Answer 2

[Jon]

I wasn't completely clear from your description what output you wanted, but I think this is what you were looking for. Still loops through rows of M, but maybe more efficient that what you had tried. Would have to time it to see.

% Build example input matrix
M = zeros (5,8) -1;
M(1,1:3) = 1:3;
M(3,2:3)  = 2:3;
M(3,5:7) = 5:7;
M(4,2:5) = 2:5;
M

M = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

% Make matrix that has ones where there are values and 0 otherwise
V = M ~= -1;

% Loop through rows finding midpoint values of longest blocks
[m,n] = size(V);
midVals = NaN(m,1); % preallocate array to hold middle region values
for k = 1:m
    % Find indices of beginning and end of each block of values
    idx = [true,diff(V(k,:))~=0,true];

    % Find length of each block
    blklngth = diff(find(idx));

    % Make a vector whose elements give length of the block that each
    % element in original matrix belongs to
    x = repelem(blklngth,blklngth).*V(k,:);

    % Find the index of the middle (rounding up) of the longest block
    % and assign middle value to output vector
    maxBlkLngth = max(x);
    if maxBlkLngth > 0
        idxMid = round(mean(find(x == max(x))));
        midVals(k) = M(k,idxMid);
    end

end
midVals

midVals = 5×1

2 NaN 6 4 NaN

[Jon]

I should also credit @Jan for his posting in https://www.mathworks.com/matlabcentral/answers/382011-how-to-count-the-number-of-consecutive-identical-elements-in-both-the-directions-in-a-binary-vecto from which I borrowed a lot of ideas used in the above

[Jon]

I've cleaned up my looping through rows approach, which I will first demonstrate for your simple example, and then I compare timing on a big array with Bruno's no loop approach.

% Build example input matrix
M = zeros (5,8) -1;
M(1,1:3) = 1:3;
M(3,2:3)  = 2:3;
M(3,5:7) = 5:7;
M(4,2:5) = 2:5;
M

M = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

% Make matrix that has ones where there are values and 0 otherwise
V = M ~= -1;

% Loop through rows finding midpoint values of longest blocks
[m,n] = size(V);
midVals = NaN(m,1); % preallocate array to hold middle region values
for k = 1:m
    % Find indices of start and end of each block of values, if there are
    % any
    if any(V(k,:))
        idx = find(diff([false,V(k,:),false]));
        iStrt = idx(1:2:end);
        iEnd = idx(2:2:end);

        % Find number of  elements in each block
        nBlk = iEnd - iStrt;

        % Find the largest block,and its position in the
        % list of blocks
        [nMax,iPos] = max(nBlk);

        % Find the midpoint (rounding up) index of the longest block
        idxMid = iStrt(iPos) + round((nMax-1)/2);

        % Assign the midpoint value
        midVals(k) = M(k,idxMid);
    end
end
midVals

midVals = 5×1

2 NaN 6 4 NaN

Now compare timing on @Bruno Luong's big matrix, using my loop through rows, and Bruno's no loop method (Note, I provide the value of the midpoint element for each row. I'm unclear about the meaning of first_idx, second_idx in Bruno's approach, but I left it as is, hopefully this is still a fair comparison)

M = repmat(1:1000, 10000, 1);
M(randi(end, 1000000, 1)) = -1;

tic

% Make matrix that has ones where there are values and 0 otherwise
V = M ~= -1;

% Loop through rows finding midpoint values of longest blocks
[m,n] = size(V);
midVals = NaN(m,1); % preallocate array to hold middle region values
for k = 1:m
    % Find indices of start and end of each block of values, if there are
    % any
    if any(V(k,:))
        idx = find(diff([false,V(k,:),false]));
        iStrt = idx(1:2:end);
        iEnd = idx(2:2:end);

        % Find number of  elements in each block
        nBlk = iEnd - iStrt;

        % Find the largest block,and its position in the
        % list of blocks
        [nMax,iPos] = max(nBlk);

        % Find the midpoint (rounding up) index of the longest block
        idxMid = iStrt(iPos) + round((nMax-1)/2);

        % Assign the midpoint value
        midVals(k) = M(k,idxMid);
    end
end
toc

Elapsed time is 0.163425 seconds.

% Bruno's no loop method
tic
A = (M>0).';
n = size(A,2);
f = false(1, n);
[i,j] = find(diff([f; A; f]));
c = i(1:2:end);
w = i(2:2:end)-c;
r = j(1:2:end);
[rw,is] = sortrows([r,w],[1 -2]);
keep = [true; diff(rw(:,1))>0];
rw = rw(keep,:);
c = c(is(keep));
r = rw(:,1);
w = rw(:,2);
first_idx = accumarray(r, c-1+ceil(w/2), [n,1])';
second_idx = logical(first_idx);

toc

Elapsed time is 0.261603 seconds.

So it seems like the looping through rows is not that expensive, actually time is less (maybe I am missing something?) And while, there is no doubt that Bruno's use of indexing to achieve a no loop approach is very elegant, at least for me the row by row approach is more obvious.

[Dominik Rhiem]

Thank you a lot for your work, I appreciate it. I might test it out, but there is a simple reason for me to not want to rely on a loop here (which I probably also should have mentioned in the beginning): I have a GPU available which I would like to use to speed up calculations, and in my experience, a GPU is really slow with a loop. I have tested Bruno's approach with and without GPU support, and in my test cases, the GPU speeds up calculations by a factor of 3-4 compared to pure CPU, while your example reaches a factor of 2 at best while being practically restricted to CPU (don't get me wrong, this is really impressive to me!). I agree, though, that the row by row approach is more obvious.

[Jon]

I haven't worked with GPU's so I can't offer any guidance on that. Surprising though that for a simple loop, as in my code, where each iteration is independent of previous iterations, that the iterations of the loop can't be split over multiple processing units and run in parallel. I have read a little about parfor, which I guess does this on multiple cpu cores. Not sure if there might be some equivalent for GPU's.

[Bruno Luong]

The loop is surprisingly fast. I know it but I'm still puzzeling each time.

[Jon]

In case you're still interested, here's another no-loop approch. The idea is to unwrap the matrix columwise and find all of the block run lengths in one go. I think the approach is relatively straightforward and easy to understand.

Benchmarking this, it seems that my new no-loop approach runs slightly faster than either my previous loop approach, or Bruno's no-loop approach. You would have to see whether this holds up when running on GPU's. Note however that since the example data is generated randomly, the benchmark times can vary. I also noticed that there was more variability in the durations, and in the ordering of the durations (which one was the fastest) for the three approaches running it here on the MATLAB answers server, and on my own computer. In any case they are all roughly the same.

% Build example input matrix
M = zeros (5,8) -1;
M(1,1:3) = 1:3;
M(3,2:3)  = 2:3;
M(3,5:7) = 5:7;
M(4,2:5) = 2:5;
M

M = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

% Get dimensions for later use
[m,n] = size(M);

% Make matrix that has ones where there are values and 0 otherwise
V = (M > 0);
V = [V false(m,1)]'; % Provide break in blocks at end of rows

% Viewing matrix values as one long column, find beginning and end of each
% block of values
idx = find(diff([false;V(:);false]));

% Find number of  elements in each block
iStrt = idx(1:2:end);
iEnd = idx(2:2:end);
nBlk = iEnd - iStrt;

% Store the block length in the starting element of each block
% (can't store length values in V because it is a logical matrix)
L = zeros(size(V));
L(iStrt) = nBlk;

% Viewing L  as a matrix, find longest block in each column, and its
% location
[nMax,iPos] = max(L);

% Find the midpoint (rounding up) row index of the longest block
idxMid = iPos + round((nMax-1)/2);

% Only use columns (rows in original M) that hve blocks of values
idl = idxMid ~= 0;
cols = find(idl);

% Assign the midpoint value, for rows in M that have no blocks, assign NaN
midVals = NaN(m,1);
idx = sub2ind([m,n],cols,idxMid(idl));  % Linear indices to get values
midVals(cols) = M(idx);


midVals

midVals = 5×1

2 NaN 6 4 NaN

Here are the benchmark runs:

M = repmat(1:1000, 10000, 1);
M(randi(end, 1000000, 1)) = -1;


% My no-loop method

tic
% Get dimensions for later use
[m,n] = size(M);

% Make matrix that has ones where there are values and 0 otherwise
V = (M > 0);
V = [V false(m,1)]'; % Provide break in blocks at end of rows

% Viewing matrix values as one long column, find beginning and end of each
% block of values
idx = find(diff([false;V(:);false]));

% Find number of  elements in each block
iStrt = idx(1:2:end);
iEnd = idx(2:2:end);
nBlk = iEnd - iStrt;

% Store the block length in the starting element of each block
% (can't store length values in V because it is a logical matrix)
L = zeros(size(V));
L(iStrt) = nBlk;

% Viewing L  as a matrix, find longest block in each column, and its
% location
[nMax,iPos] = max(L);

% Find the midpoint (rounding up) row index of the longest block
idxMid = iPos + round((nMax-1)/2);

% Only use columns (rows in original M) that hve blocks of values
idl = idxMid ~= 0;
cols = find(idl);

% Assign the midpoint value, for rows in M that have no blocks, assign NaN
midVals = NaN(m,1);
idx = sub2ind([m,n],cols,idxMid(idl));  % Linear indices to get values
midVals(cols) = M(idx);

toc

Elapsed time is 0.184442 seconds.

% My Earlier Loop Method
tic
% Make matrix that has ones where there are values and 0 otherwise
V = M ~= -1;

% Loop through rows finding midpoint values of longest blocks
[m,~] = size(V);
midVals = NaN(m,1); % preallocate array to hold middle region values
for k = 1:m
    % Find indices of start and end of each block of values, if there are
    % any
    if any(V(k,:))
        idx = find(diff([false,V(k,:),false]));
        iStrt = idx(1:2:end);
        iEnd = idx(2:2:end);

        % Find number of  elements in each block
        nBlk = iEnd - iStrt;

        % Find the largest block,and its position in the
        % list of blocks
        [nMax,iPos] = max(nBlk);

        % Find the midpoint (rounding up) index of the longest block
        idxMid = iStrt(iPos) + round((nMax-1)/2);

        % Assign the midpoint value
        midVals(k) = M(k,idxMid);
    end
end
toc

Elapsed time is 0.222911 seconds.

% Bruno's no-loop method
tic
A = (M>0).';
n = size(A,2);
f = false(1, n);
[i,j] = find(diff([f; A; f]));
c = i(1:2:end);
w = i(2:2:end)-c;
r = j(1:2:end);
[rw,is] = sortrows([r,w],[1 -2]);
keep = [true; diff(rw(:,1))>0];
rw = rw(keep,:);
c = c(is(keep));
r = rw(:,1);
w = rw(:,2);
first_idx = accumarray(r, c-1+ceil(w/2), [n,1])';
second_idx = logical(first_idx);

toc

Elapsed time is 0.244201 seconds.

[Bruno Luong]

Not bad, the new neat idea to me is push back

L = zeros(size(V));
L(iStrt) = nBlk;

[Jon]

Thanks. Wish I could store the blockLengths (nBlk in above) in the V matrix instead of having to create a whole new one. Trouble is that since V is a logical matrix, it can't store the values. I could convert the V matrix to a numerical matrix (double or uint##) but not sure whether this would be any better (I guess it would use a little less space, but I'm not sure about performance).

[Bruno Luong]

@Jon I modify your code so that you don't have to allocate another big array

% Get dimensions for later use
[m,n] = size(M);

% Make matrix that has ones where there are values and 0 otherwise
V = (M > 0);
V = [V false(m,1)]'; % Provide break in blocks at end of rows

% Viewing matrix values as one long column, find beginning and end of each
% block of values
V = diff([false;V(:)]); % Bruno's comment: no need last false since V altrady has false at the end
                        % Another reason, it preserve V to have same number
                        % of elements, that allow me to reshape later
idx = find(V);

% Find number of  elements in each block
iStrt = idx(1:2:end);
iEnd = idx(2:2:end);
nBlk = iEnd - iStrt;

% Store the block length in the starting element of each block
% (can't store length values in V because it is a logical matrix)
V = reshape(V,[n+1,m]); % Bruno's comment: This was L, I reuse V which has 1 at iStrt 'replaced) and -1 at iEnd
V(iStrt) = nBlk;

% Viewing L  as a matrix, find longest block in each column, and its
% location
[nMax,iPos] = max(V);

% Find the midpoint (rounding up) row index of the longest block
idxMid = iPos + round((nMax-1)/2);

% Only use columns (rows in original M) that hve blocks of values
idl = idxMid ~= 0;
cols = find(idl);

% Assign the midpoint value, for rows in M that have no blocks, assign NaN
midVals = NaN(m,1);
idx = sub2ind([m,n],cols,idxMid(idl));  % Linear indices to get values
midVals(cols) = M(idx);

[Jon]

Hi Bruno, thanks for your suggested modification. The problem is (and why I needed to make the numerical, matrix L, rather than storing in the logical matrix V) that assigning a non-zero double, to an element of a logical matrix, just sets the value to 1 (true). So for example

V = false(2,3)

V = 2×3 logical array
   0   0   0
   0   0   0

V(2,2) = 23

V = 2×3 logical array
   0   0   0
   0   1   0

Therefore, I don't think the reshaping you do for V solves the problem. I could already assign to the correct element of V using the linear indices, even without reshaping. Problem is value gets stored as a 1

[Bruno Luong]

But I change V to double just before find (you don't need the original V afterwards).

V metamorphoses and has three lifes

[Jon]

Sorry, Bruno, I didn't read your modified file carefully enough. Also missed that by reassigning V to diff([false;V(:)]) V has now become a double vector. So now I see that as a double we can assign the block lengths to it, and because it is a vector we must later reshape it. That's great! Thanks so much. Your code is always inspiring.

Answer 3

[Mathieu NOE]

hello

I tried to put some code together and ended with that ; maybe interesting (?)

it will detect the largest region for each row and store (and display - see red crosses ) the center position of each region

I haven't rounded these values so add it if you need it

matrix = zeros (5,8) -1;
matrix(1,1:3) = 1:3;
matrix(3,2:3)  = 2:3;
matrix(3,5:7) = 5:7;
matrix(4,2:5) = 2:5;
matrix

matrix = 5×8

1 2 3 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 3 -1 5 6 7 -1 -1 2 3 4 5 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

%% main code 
mat = (matrix>0); % creat logical array 

for k = 1:5
    
    [begin,ends] = find_start_end_group(mat(k,:));
    ll = ends-begin+1;
    if isempty(ll)
        index_pos(k) = NaN;
    else
        if numel(ll)>1 % more than one region in this row
            [v,id] = max(ll);
            index_pos(k) = 0.5*(begin(id)+ends(id)); % % center position (index) of largest region
        else
            index_pos(k) = 0.5*(begin+ends); % % center position (index) of largest region
        end
    end

end

% plot
cmap = [0 0 0; parula(128)];
matrix(matrix<0) = NaN;
imagesc(matrix)
colormap(cmap);
caxis([-1 10]);
colorbar('vert');
hold on
for k = 1:5
    plot(index_pos(k),k,'r+','markersize',25);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [begin,ends] = find_start_end_group(ind)
% This locates the beginning /ending points of data groups
% Important : ind must be a LOGICAL array

    D = diff([0;ind(:);0]);
    begin = find(D == 1);
    ends = find(D == -1) - 1;
end

Answer 4

[Fifteen12]

I've also put some work into this, but I'm unclear with what you want as your output, so I've got it to a good stopping place and I'll let you finish it with what you want to do. This script basically identifies the first and last index for every region in an inputted matrix using vectorized calls.

function [first, last] = findRegion(in)
    valid = find((in + 1)');
    neighbors = [0; valid];
    regions = [valid; numel(in) + 1] - [-1; valid]; %Spaces between each valid index (first cell inflated for easier processing, last cell is handled later)
    starting_indices = valid(regions > 1);
    temp = logical([regions > 1; 0]);
    ending_indices = valid(temp(2:end));
    if (starting_indices(end)) == numel(in)
        ending_indices(end+1) =  numel(in);
    end
    
    %% Separate regions on different rows
    rows = transpose(0:width(in):numel(in));
    row_pairs = [floor((starting_indices - 1) ./ width(in)), floor((ending_indices - 1) ./ width(in))];
    wrapped_regions = logical(row_pairs(:, 1) - row_pairs(:, 2));
    first = sort([starting_indices; rows(wrapped_regions) + 1]);
    last = sort([ending_indices; rows(wrapped_regions)]);
end

I wasn't sure if the largest region was the one with the most values or the one with the highest aggregated value, but I think you can probably calculate either from here.

Answer 5

[Image Analyst]

If you have the Image Processing Toolbox, it's pretty easy:

M = [...]
    1     2     3    -1    -1    -1    -1    -1
    -1    -1    -1    -1    -1    -1    -1    -1
    -1     2     3    -1     5     6     7    -1
    -1     2     3     4     5    -1    -1    -1
    -1    -1    -1    -1    -1    -1    -1    -1]
[rows, columns] = size(M); % Get dimensions.
% Process each row one at a time.
for row = 1 : rows
    mask = M(row, :) ~= -1;
    props{row} = regionprops(mask, 'PixelIdxList')
    % Now, not sure how to define the "central value rounded up"
end

So I'm just not sure how you define central value. What would be the central value of a vector like this: [3,6,1,9]?

[Image Analyst]

I know you were happy with your times of around 0.2 seconds, but with my solution I'm getting around 25 times faster than that on my 12 year old computer:

% Process each row one at a time.
tic
for row = 1 : rows
    mask = M(row, :) ~= -1;
    props{row} = regionprops(mask, 'PixelIdxList');
    % Now, not sure how to define the "central value rounded up"
end
toc

Elapsed time is 0.007568 seconds.

However it does require the Image Processing Toolbox for the regionprops function, which is specifically made to detect and measure regions like that.

[Bruno Luong]

@Image Analyst 0.1-0.2 second is time with data of size 1000 x 10000 as bellow, not with OP toy example. Your code take 4 seconds which is more than 20 time slower than our codes (without tollbox).

M = repmat(1:1000, 10000, 1);
M(randi(end, 1000000, 1)) = -1;

tic
for row = 1 : size(M,1)
    mask = M(row, :) ~= -1;
    props{row} = regionprops(mask, 'PixelIdxList');
    % Now, not sure how to define the "central value rounded up"
end
toc

Elapsed time is 4.032186 seconds.