2

I've got an image in matlab defined on a real scale (0,1) and a mask defined on a integer scale.

Example

mask = [ 1 1 1 3 4 ;
         1 1 1 2 4 ;
         1 1 2 2 2  ]

img = [ 0.1 0.1 0.2 0.2 0.3 ;
        0.1 0.1 0.2 0.3 0.3 ;
        0.1 0.1 0.3 0.3 0.3 ]

and for each region in the mask (i.e. 1,2,3,4) I want to compute a certain feature (say, the mean) on the corresponding image's intensities.

The algorithm I've used is 

for i = labels
  region = img(mask==i);
  feature(i) = mean(region);
end

Now, this algorithm is pretty slow for images of size 300x400x500, and a cardinality of the labels set > 40000 (which, btw, is exactly my case).

Any suggestion on how to speed up my code?

user1384636
  • 481
  • 6
  • 17
  • What do you mean by *images of size `300x400x500`*? Is it 500 images of size `300x400`? – knedlsepp Jan 28 '15 at 15:04
  • 1
    @knedlsepp It is obviously a volumetric image, a 3-dimensional array in Matlab. – Jonathan H Jan 28 '15 at 15:06
  • 2
    @Sh3ljohn: I am asking because the example is not in line with the data given. The OP could either have: 2D `mask`, 3D `img` and wants the same mask for all the layers of `img` along the third dimension, or have a 3D `mask` and 3D `img`. That's why it isn't *that* obvious. – knedlsepp Jan 28 '15 at 15:09
  • @knedlsepp The lack of relevant information in the OP is shocking given the question asked. But to summarize the comments, the OP is doing clustering in a small volumetric image, and wishes to compute average, sdev and skewness of clusters intensities as fast as possible in Matlab. – Jonathan H Jan 28 '15 at 15:17
  • @Sh3ljohn: Oh! In the comment to your answer! So `mask` really is a 3D array. – knedlsepp Jan 28 '15 at 15:21
  • sorry if I didn't give an example with volumetric data, I thought It was pretty clear what I was asking – user1384636 Jan 28 '15 at 18:14
  • For someone in the image processing field, it is obvious, but for those who may not be, it isn't. When posting questions in the future.... just a suggestion... assume that everyone is an idiot. Give as much detail as you can, with everything you have tried to solve your problem. Code that actually runs will give you brownie points. – rayryeng Jan 28 '15 at 19:41

6 Answers6

3

The regionprops function in the Image Processing Toolbox should do it for you. For example, use this syntax:

stats = regionprops(L,I, 'MeanIntensity');

to get the mean value for every region. L is the array containing the labels. I is your image.

Dima
  • 38,860
  • 14
  • 75
  • 115
1

bsxfun and fast matrix-multiplication based approach to get the mean values -

num_labels = max(mask(:)); %// number of labels used
labels = 1:num_labels;     %// array of all labels

%// Get a 2D mask for all labels, with each column representing a label.
%// This is a setup for use with matix-multiplication later on
mask_labels = bsxfun(@eq,mask(:),labels);

%// Perform fast matrix multiplication as a way to perform summation for
%// all labels and then do elementwise division to get the mean values
feature_out = (img(:)'*mask_labels)./sum(mask_labels,1);
Community
  • 1
  • 1
Divakar
  • 218,885
  • 19
  • 262
  • 358
1

To get the mean elegantly (though it may not be the absolute fastest approach), you can use accumarray:

meanFeatures = accumarray(mask(:),image(:),[],@mean);

If you want e.g. mean and std, you use

meanAndStd = cell2mat(accumarray(mask(:),image(:),[],@(x){mean(x), std(x)}));
Jonas
  • 74,690
  • 10
  • 137
  • 177
  • I think going for `accumarray` is the logical choice (with the possible addition of `unique` to get safe masks with entries in `1:length(unique(mask(:))`), never mind iterating over the images using a loop. – knedlsepp Jan 28 '15 at 15:17
  • Oh! `mask` is supposed to be a 3D array. Never mind. – knedlsepp Jan 28 '15 at 15:22
0

I think we need to know more about:

  • How you compute the mask?
  • What are the features that you want to compute on these regions? (Mean only?)
  • What do you mean by slow?

If the formula that assigns values to each element of the mask can be vectorized, then so can the access to your volume.

In the general case though, I would say that this would be easy to do using Mex files instead. If you want to go down that way, I could recommend using my tiny library to ease the transition to Mex, but you'll need to compile using C++11 (-std=c++0x flag).

If you want to stay in Matlab, and that you only want to compute averages on these regions, then accumarray is the function you are looking for.

Jonathan H
  • 7,591
  • 5
  • 47
  • 80
  • How you compute the mask? It's a procedure to create supervoxels, but is actually not important how I create it, let's say that it has been given. What are the features that you want to compute on these regions? (Mean only?) Mean, standard deviation of intensities, skewness. What do you mean by slow? > 10 minutes. I repeat that I have an image of size 300x400x500 and the corresponding mask has (of course) the same size. – user1384636 Jan 28 '15 at 14:25
  • I asked because it _is_ important. Without further details about the procedure that creates the mask, and given the set of features you wish to compute, Mex files is the way to go. That said, given the simplicity of the features you compute and the small size of your volume, your Matlab code is probably poorly written; then again, without more information about this code, I can't say why it is so slow. – Jonathan H Jan 28 '15 at 14:52
  • @user1384636: Your comment here should really be added to the question itself! – knedlsepp Jan 28 '15 at 15:24
0

To the best of my knowledge, this is not the kind of problem Matlab is very good at. I would try the naive - non Matlab approach:

 label_sum = zeros(1,num_labels);
 label_size = zeros(1,num_labels);
 for z=1:500,
    for x=1:400,
       for y=1:300,
          label = labels(y, x, z);
          label_sum(label) += img(y, x, z);
          label_size(label) += 1;
        end
     end
  end
  label_sum(label_size > 0) ./ label_size(label_size > 0)

it is not ideal but will probably take way less than 10 minutes...

Rosa Gronchi
  • 1,828
  • 15
  • 25
  • maybe, but it is a very simple one to try and Matlab has become much better in performing loops in recent years. – Rosa Gronchi Jan 28 '15 at 14:56
  • It does not answer the question though, unless the OP implementation is very poorly written, your code is pretty much the slowest solution I could think of, and certainly slower than the algorithm given in the OP (again, if written correctly). – Jonathan H Jan 28 '15 at 14:57
  • in this implementation I access each element in the 3D image once. The implementation suggested in the question access each element in the 3D image once per label (== 40,000 times), so I think it worth trying. – Rosa Gronchi Jan 28 '15 at 15:02
  • 1
    The _algorithm_ you use is efficient; in fact, the most efficient Mex/C++ solution I can think of is very similar to this. But Matlab won't execute that code efficiently. – Jonathan H Jan 28 '15 at 15:03
0

You may use arrayfun and avoid using the for loop. Using your labels vector of values to be iterated, it goes as:

features=arrayfun(@(x) mean(img(mask==x)), labels);

The syntax @(x) will iterate upon each element of labels, replacing in this way the for loop. You may also use any function instead of mean, either it is a built-in function or a function that you have created. You should check its speed, but I think this is the fastest way.

fyts
  • 493
  • 5
  • 12
  • One would hope so. Yet in reality, `arrayfun` is even slower than the original `for`-approach... It may be a nifty one-liner, but using `arrayfun` won't be a performance improvement. – knedlsepp Jan 28 '15 at 15:39