Error while finding mean and variance using blockproc of an image

I have a blockproc function as follows:

function [mean_ch,variance] = block_process(subchannel)
fun1 = @(block_struct) mean(subchannel(:));
mean_ch = blockproc(subchannel, [8 8], fun1);
fun2 = @(block_struct) var(subchannel(:));
variance = blockproc(subchannel, [8 8], fun2);
end

When I run this for, say the red channel of my image(in tiff format), I get a matrix of all same values. When I try to do a scatterplot by reshaping these values I get a single value. I believe the function is only running over a single block hence the single value, but I am supposed to use blockproc as a sliding window operator to calculate the sample mean and variance of each 8x8 block. How can achieve the right results?

Answer 1

[Walter Roberson]

function [mean_ch,variance] = block_process(subchannel)
    COL = @(M) M(:);
    fun1 = @(block_struct) mean(COL(block_struct.data(:,:,subchannel)));
    mean_ch = blockproc(subchannel, [8 8], fun1);
    fun2 = @(block_struct) var(COL(block_struct.data(:,:,subchannel)));
    variance = blockproc(subchannel, [8 8], fun2);
end

This assumes that you want to process each independent 8 x 8 block, rather than wanting do sliding window.

If you do want to do sliding window, and you want it to proceed one pixel at a time, then there is a bit of a problem because of the way that blockproc handles sliding windows: when you want to slide by an odd amount (such as 1 pixel at a time) then it can only handle blocks that are odd-sized (such as 7 x 7, but not 8 x 8)

Because of this, if you do want to shift one pixel at a time and you want 8 x 8 blocks, it is usually easier to do some conv2() operations, or nlfilter()

[Deepika Sundresh]

The problem says "Analyse each subchannel of the image with a sliding window operator that outputs the local sample means and local sample variances for each position of the window". Does this mean I cannot use blockproc?

[Walter Roberson]

You have a couple of different possibilities here:

1. You could slide by more than one pixel at a time while using an even block size such as 8 x 8; OR
2. You could slide by one pixel at a time while using an odd block size such as 7 x 7 with blockproc; OR
3. You could use nlfilter() or some careful calculations with conv2() to calculate the values you want

[Deepika Sundresh]

I made this function as a sliding window operator to compute mean and variance of the 8x8 blocks of the channel.

function [m, v] = window(channel, h, w)
  [r,c] = size(channel);

  for i = 1:c-h
      for j = 1:r-w
          M = channel(i:(i+h)-1,j:(j+w)-1);
          m(i) = mean(M(:));
          v(i) = sum(double(M(:))-repmat((m(i)),h*w,1))/(h*w);
      end
  end
  figure,
  scatter(m,v);
  end

When I run this for the channel (r = 1380, c = 1842, h and w = 8) I get the error

  Index in position 1 exceeds array bounds (must not exceed 1380)

Any idea what could b causing the problem?

[Walter Roberson]

  [r,c] = size(channel);

r reflects rows, c reflects (columns times number of image planes)

   for i = 1:c-h

i involves c, so i must be intended to index columns -- the second index.

      for j = 1:r-w

j involves r, so j must be intended to index rows -- the first index.

        M = channel(i:(i+h)-1,j:(j+w)-1);

but i is used for the first index, and j is used for the second index.

Perhaps you have incorrectly defined i and j

[Deepika Sundresh]

Ah yes, that was the reason. Thank you. Can you guide me how I can check if the scatter plot of a channel is correct or not? I obtained the scatterplot from this function but I'm unable to verify it.

[Walter Roberson]

Question: for your purposes, is it acceptable that your output is smaller than your input? Because there are ways to code it so that it operates on whatever is inside the sliding window -- though definitions of what is inside the window are certainly easier when the window is an odd size instead of an even size.

[Deepika Sundresh]

The output would be smaller than the input right? That's what I understood :)

Answer 2

[Image Analyst]

To get the mean in a sliding window that slides over one pixel at a time:

windowWidth = 3; % Whatever...
kernel = ones(windowWidth) / windowWidth^2;
blurredImage = imfilter(grayImage, kernel);

To get the standard deviation, use stdfilt():

windowWidth = 3; % Whatever...
kernel = ones(windowWidth);
sdImage = stdfilt(grayImage, kernel);
% Square to get variance
varImage = double(sdImage) .^ 2;