Most efficient way to get consequent IDs in a large matrix

Question

After performing a watershed of a 3D CT, I choose only the particles which are smaller than a certain value and larger than another one. As a final output I however need a matrix only containing consecutive IDs. My implementation is the following:

% Get unique IDs
grain_id = unique(L);

% Get rid of artefacts
% Compute histogram for each ID
% and compare volume numv with thresholds

% Reject grains smaller or larger than threshold
reject = grain_id(or(numv<vmin,numv>vmax));

% Keep 0s (boundaries) and 1 (voids)
reject = reject(3:end);        

% Rejected become void
L(ismember(L,reject))=1;

% Get number of grains
grain_id = unique(L);
numgrains = numel(grain_id);

% Consecutive IDs
idx = false(size(L));
for i=1:numel(reject)
    idx = L>reject(numel(reject)+1-i);
    L = L-uint16(idx);
end

I have a 1226x1226x3600 matrix, so performance really matters. One loop takes approx. 5s. It is very likely not the most efficient way to achieve the goal, but at the moment I do not have better ideas. Do you?

Answer 1

If I am understanding your problem statement correctly, this is precisely what the third output of the unique function is for. It examines your data and for each element, the third output provides a mapping that tells you which index into the output of unique you would need for the corresponding element. Coincidentally, this provides a new integer mapping that is consecutive from 1 up to as many unique labels as there are in the input.

However, this is returned as a column unwrapped vector for anything other than a vector as the input into unique so you'll need to reshape this back into the same dimensions as the input at the end.

Therefore, what you would need is simply:

[~,~,id] = unique(L);
id = reshape(id, size(L));

id would be the matrix that is the same size as L that you used as input into L.

Here's a toy example to be sure we are on the same page:

>> rng(123); L = randi(50, 10, 10)

L =

    35    18    32     5    32     7    34     5    16    36
    15    37    43    22     6    42    30    45    35    50
    12    22    37    22    16    31    32    32    28    18
    28     3    31    25    21    28    34    37    20    39
    36    20    37    22    44    18    43     1    47    30
    22    37    17    16    13    16     5    30    43    35
    50    10    19    22    25    21    39    28    18     8
    35     9    12    45    50    35    13     8     3    20
    25    27    15    48    26    44    10     8    16    13
    20    27    32    26    31    26    29    35    20    18

Here I create a 10 x 10 matrix of random numbers that are not consecutive. We can see this by viewing all of the unique numbers in this matrix:

>> unique(L).'

ans =

  Columns 1 through 19

     1     3     5     6     7     8     9    10    12    13    15    16    17    18    19    20    21    22    25

  Columns 20 through 38

    26    27    28    29    30    31    32    34    35    36    37    39    42    43    44    45    47    48    50

Notice the jump from 1 to 3, or from 13 to 15 for example. By performing the code I wrote above, we now get:

>> id

id =

    28    14    26     3    26     5    27     3    12    29
    11    30    33    18     4    32    24    35    28    38
     9    18    30    18    12    25    26    26    22    14
    22     2    25    19    17    22    27    30    16    31
    29    16    30    18    34    14    33     1    36    24
    18    30    13    12    10    12     3    24    33    28
    38     8    15    18    19    17    31    22    14     6
    28     7     9    35    38    28    10     6     2    16
    19    21    11    37    20    34     8     6    12    10
    16    21    26    20    25    20    23    28    16    14

As you can see here, the label 3 becomes label 2 to ensure the labelling is consecutive. Similarly, labels 13 and 15 become 10 and 11 to ensure the consecutive order you require and so on. To be absolutely sure, here is a list of all of the unique values in the output:

>> unique(id).'

ans =

  Columns 1 through 19

     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19

  Columns 20 through 38

    20    21    22    23    24    25    26    27    28    29    30    31    32    33    34    35    36    37    38