How to 'walk' a mask/ROI over an image?

Hi,

At the moment, I am working at a script to determine the homogeneity of an image. For this, I want to 'move' an rectangular ROI (for example, 0.5cm x 0.5cm) over the image in steps of half ROI. I would like to determine for example, the mean pixel value for each ROI-selection.

So far, I can not really find a useful feature or script for this. I was wondering how this can be done the best and fastest?

Thanks in advance for the reply.

[Adam]

doc nlfilter

should be able to do this, although it does come with an embedded progress bar which is not always desirable.

Answer 1

[Walter Roberson]

blockproc(), and make sure you specify BorderSize. Be careful with TrimBorder when you use overlaps.

[Kevin Semmelink]

This function works, but I don't see how I am able to force the block to move with steps of half the size of the ROI. I've found an example online to specify what I want.

[Walter Roberson]

H = 2; W = 2; OH = 2; OW = 2;
output = blockproc(YourImage, [H W], @(block) ones(size(block.data,1),size(block.data,2))*mean2(block.data), 'BorderSize', [OH/2 OW/2], 'TrimBorder', 1)

This would pass in to the anonymous function a data block that is

((OH/2) + H + (OH/2)) by ((OW/2) + W + (OW/2))

which would be 1+2+1 by 1+2+1 = 4 x 4. Then after the calculation for that 4 x 4 area, the focus of the block would be moved H or W further along.

The action done for each 4 x 4 block would be to take the mean2() of the block, and replicate the result to the same size as the block. Then, because TrimBorder is set, blockproc would trim off OH/2 from the top and OH/2 from the bottom and OW/2 from the left and OW/2 from the right and store the resulting 2 x 2 in the array.

Note: watch out for partial blocks; see the PadPartialBlocks setting. The code I use here automatically chooses the correct output size and the TrimBorder setting knows how to deal with trimming the partial blocks correctly, so the behavior for what I show should be okay.

The function can be simplified slightly if you know your image is 2D instead of 3D.

The output will have 2 x 2 blocks of the same value. That is because you did not define the size of output you want. mean of a 2 x 2 block is a scalar but it is not clear that you wanted a scalar result per block as the output for that would be half the size of the input in rows and columns. The equivalent call that gives the smaller output would be

H = 2; W = 2; OH = 2; OW = 2;
output = blockproc(YourImage, [H W], @(block) mean2(block.data), 'BorderSize', [OH/2 OW/2], 'TrimBorder', 0)

You could also consider a different approach:

 np = size(YourImage,3);
 output = convn(YourImage, ones(4,4,np)./(4*4*np), 'same');
 output = output(1:2:end, 1:2:end, 1);

What this does is walk a 4 x 4 x number_of_panes mask over the image one pixel at a time, and then throw away half of the results in rows and columns and keep only one plane.

If you know for sure that your image is 2D then you can simplify to

output = conv2(YourImage, ones(4,4)./16, 'same');
output = output(1:2:end, 1:2:end);

[Walter Roberson]

Note:

To move "half" an ROI over each time, then it follows immediately that the ROI length must even so that half of the length is integer.

But then, because the remaining half of the ROI is divided between the two sides, the remaining half of the ROI must also be even.

Therefore to use blockproc for this purpose, it is necessary that the original ROI length must be a multiple of 4.

The restriction that the ROI length must be a multiple of 4 does not hold for using convolution to doing the calculation. It does still need to be even, though, for "move by half" to have meaning.

For the 2D version:

RL = 10;     %full ROI length
half = RL/2;

output = conv2(YourImage, ones(RL,RL)./(RL*RL), 'same');
output = output(1 : half : end, 1 : half : end);

Note: watch out for edge effects on the convolution.

[Kevin Semmelink]

The last step you wrote works like a charm! I am not yet sure why, but it worked after I changed 'same' to 'valid' as shown below. After that, it determines exactly the mean pixel value of each 'half' ROI-step.

 output = conv2(YourImage, ones(RL,RL)./(RL*RL), 'valid');
 output = output(1 : half : end, 1 : half : end);

I tested it with a random 16x16 matrix and calculated the mean pixel values step by step and with use of the convolution. This resulted in two 7x7 matrix with corresponding values!

Unfortunately, the code with the blockproc did not correspond. Do you have any idea why? I would like to determine things like the standard deviation next. I think, this is easier when I could use the function blockproc, since I only have to change mean2 to std2. Am I right?

[Kevin Semmelink]

The manual test was something like this:

 test_1 = A(1:4,1:4);
 test_2 = A(1:4,3:6);
 test_3 = A(1:4,5:8);
 test_4 = A(3:6,1:4);
 test_5 = A(3:6,3:6);
 test_6 = A(3:6,5:8);

 ave1 = mean2(test_1(:));
 ave2 = mean2(test_2(:));
 ave3 = mean2(test_3(:));
 ave4 = mean2(test_4(:));
 ave5 = mean2(test_5(:));
 ave6 = mean2(test_6(:));

 FIN = [ave1 ave2 ave3; ave4 ave5 ave6];

This 2x3 matrix was corresponding with

output(1:2,1:3)

A = my image

 A = rand(16,16);

[Walter Roberson]

'valid' would be a step towards dealing with the edge effects I mentioned.

It is easier to use blockproc for standard deviation as you have noted, but as I explained above if you want your ROI to advance by half the width then your ROI would have to be a multiple of 4 in length, which would not allow an ROI length of 10 x 10.

I just thought of a work-around that is not outright terrible (but it would be nice if there were a better way):

You can run blockproc 4 times:

p11 = blockproc(YourImage, [RL RL], AppropriateFunction );
p21 = blockproc(YourImage(RL/2+1:end, :), [RL RL], AppropriateFunction );
p12 = blockproc(YourImage(:, RL/2+1:end, [RL RL], AppropriateFunction );
p22 = blockproc(YourImage(RL/2+1:end,RL/2+1:end), [RL RL], AppropriateFunction );

That is, do four runs with no overlap, with half-window offsets.

Then, interweave the results.

%code has been tested
vert1 = reshape( [p11; p21], size(p11,1), [] );
vert2 = reshape( [p12; p22], size(p12,1), [] );

merged = reshape([vert1(:).'; vert2(:).'], 2*size(p11,1), []);

[Kevin Semmelink]

Thanks a lot! The last part didn't work on my code, but when I switched p21 and p12, it worked like a charm!

This is the final code I used.

 vert1 = reshape( [p11; p12], size(p11,1), [] );
 vert2 = reshape( [p21; p22], size(p12,1), [] );
 merged = reshape([vert1(:).'; vert2(:).'], 2*size(p11,1), []);
 merged_fin = merged(1:size(merged,1)-1,1:size(merged,2)-1);

Thanks a lot for your answers and patience!

[Walter Roberson]

You could abbreviate that last line:

 merged_fin = merged(1:end-1, 1:end-1);

[Kevin Semmelink]

You're completely right!