Thread Subject:

Difference between pixel values and centre pixel in a local neighbourhood

Hi,

My problem is as follows:

- Using nlfilter, I would like to calculate the difference between every pixel and the centre pixel of that neighbourhood (A(i) - A(x)) where A(i) is the pixel value and A(x) is the pixel value of the centre pixel
- Returning the value of the difference at each A(i)

The most trivial example for the following 2D neighbourhood:

A = magic(3)

A =

     8 1 6
     3 5 7
     4 9 2

b = floor(([3 3]+1)/2)

b =

     2 2

B = zeros(9,1);
for k = 1:numel(B)
B(k) = A(ind2sub(size(A),k)) - A(b(1),b(2));
end
B = reshape(B,3,3);
B(b(1),b(2)) = A(b(1),b(2))

B =

     3 -4 1
    -2 5 2
    -1 4 -3

Can anyone point me out how to extend this to 2 or 3D using nlfilter. I want to avoid for loops but I am not too sure how to extend this.

Thanks !

That's a linear filter, not a non-linear filter. Anyway, you don't need all that complicated stuff. All you need is conv2().

kernel = [-1 -1 -1; -1 8 -1; -1 -1 -1];
outputImage = conv2(grayImage, kernel);

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <6023a68e-9bec-4b24-adef-add7cc9d0da9@googlegroups.com>...
> That's a linear filter, not a non-linear filter. Anyway, you don't need all that complicated stuff. All you need is conv2().
>
> kernel = [-1 -1 -1; -1 8 -1; -1 -1 -1];
> outputImage = conv2(grayImage, kernel);

Hi ImageAnalyst,

Thank you for your comment.

I initially considered using a convolution, but using your specified technique does not produce the required answered.

I realise nlfilter was the wrong function to try and use, blockproc or colfilt seem more appropriate as I need to operate on distinct blocks.

I've tried the following using mat2tiles (written by Matt J) and cellfun

1) Function which calculates the difference between all pixels and the centre one
function output_nhood = diff_pix(input_nhood)
[x,y,z] = size(input_nhood);
center_pix = floor(([x y z] + 1)/2);

output_nhood = zeros(numel(input_nhood),1);
for k = 1:numel(output_nhood)
output_nhood(k) = input_nhood(ind2sub(size(input_nhood),k)) - input_nhood(center_pix(1),center_pix(2),center_pix(3));
end
output_nhood = reshape(output_nhood,x,y,z);
output_nhood(center_pix(1),center_pix(2),center_pix(3)) = input_nhood(center_pix(1),center_pix(2),center_pix(3));

2)
cellA = mat2tiles(A,n_x,n_y,n_z); % n_x, n_y and n_z are the neighbourhood dimensions
A_diff = cellfun(@(x) diff_pix(x),cellA,'UniformOutput',false);

I appreciate the fact using cells is most probably the least elegant solution, and this is represented by the fact it is taking way too long to compute, I am dealing with a large matrix.

Any pointers how to optimise this process?

Thanks !

KM

Since you're not using sliding blocks but blocks that "jump" in jumps of the window size, you'd have to use blockproc(). Be sure that your function in blockproc returns a 3x3 array, not a single number of else your output image will be 1/3 as big as your input image. If you still can't figure it out, write back.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <f4091145-322f-4184-8bb8-526f5003dbec@googlegroups.com>...
> Since you're not using sliding blocks but blocks that "jump" in jumps of the window size, you'd have to use blockproc(). Be sure that your function in blockproc returns a 3x3 array, not a single number of else your output image will be 1/3 as big as your input image. If you still can't figure it out, write back.

Hi ImageAnalyst,

I've modified my code a bit as I've realised I need to do things differently. However, I still use Blockproc in a certain section (as you've suggested) and I've been surprised by how slow it is.

Instead of using blocks from the entire image, I pre-select certain voxels based on a mask, and create blocks of size 7x7x7. The amount of voxels I pre-select is 1,602,657; thus leading to 1,602,657 (n) blocks of size 7x7x7.

I then do the following:

matrix_reshape = reshape(matrix,7^3,n);
fun = @(block_struct) diff_pix(block_struct.data);
test = blockproc(matrix_reshape,[7^3 1],fun);

where diff_pix is as previously defined in my past post, however re-written with no for loop.

Blockproc gives me an estimated time of 59 minutes ! Is there a way around this?

"K M" wrote in message <jujtk8$eat$1@newscl01ah.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <f4091145-322f-4184-8bb8-526f5003dbec@googlegroups.com>...
> > Since you're not using sliding blocks but blocks that "jump" in jumps of the window size, you'd have to use blockproc(). Be sure that your function in blockproc returns a 3x3 array, not a single number of else your output image will be 1/3 as big as your input image. If you still can't figure it out, write back.
>
> Hi ImageAnalyst,
>
> I've modified my code a bit as I've realised I need to do things differently. However, I still use Blockproc in a certain section (as you've suggested) and I've been surprised by how slow it is.
>
> Instead of using blocks from the entire image, I pre-select certain voxels based on a mask, and create blocks of size 7x7x7. The amount of voxels I pre-select is 1,602,657; thus leading to 1,602,657 (n) blocks of size 7x7x7.
>
> I then do the following:
>


> matrix_reshape = reshape(matrix,7^3,n);
> fun = @(block_struct) diff_pix(block_struct.data);
> test = blockproc(matrix_reshape,[7^3 1],fun);
>
> where diff_pix is as previously defined in my past post, however re-written with no for loop.
>
> Blockproc gives me an estimated time of 59 minutes ! Is there a way around this?

Edit: The estimated run-time fluctuates between 14-18minutes. Still, I'm sure this can be optimised ? (I dont have the Parallel Computing Toolbox)

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <f4091145-322f-4184-8bb8-526f5003dbec@googlegroups.com>...
> Since you're not using sliding blocks but blocks that "jump" in jumps of the window size, you'd have to use blockproc(). Be sure that your function in blockproc returns a 3x3 array, not a single number of else your output image will be 1/3 as big as your input image. If you still can't figure it out, write back.

I have solved the problem using some indexing and reshaping the matrix. In the end, blockproc was not needed.

Thank you for your help.

Solution:

% Random data
A = ceil(10*rand([9 12]))
% block size, note: size(A) must be multiple off w
w = [3 3];

% Engine
n = size(A);
s = [w; n./w];
B = reshape(A,s(:)');
ind = arrayfun(@(w,n) (w+1)/2:w:n, w, n, 'unif', 0);
Amid = A(ind{:});
s(1,:) = 1;
Am = reshape(Amid,s(:)');
B = reshape(bsxfun(@minus, B, Am), size(A));
B(ind{:}) = Amid;

% Check
B

% Bruno

Well you didn't use blockproc correctly. If you'd like one of my demos using it, let me know.

The code I post earlier is for block.

For sliding window, the code is as following

% Data
A=ceil(10*rand(6,8));
w = [3 3];

% Engine
n = size(A);
m = n-w+1;
dfun = @(w) linspace(-(w-1)/2,+(w-1)/2,w);
d = arrayfun(dfun, w, 'unif', 0);
D = cell(size(d));
[D{:}] = ndgrid(d{:});
B = zeros([prod(w) m]);
ind = arrayfun(@(w,m) (w-1)/2+(1:m), w, m, 'unif', 0);
Amid = A(ind{:});
dots = repmat({':'},1,n);
for p=1:numel(D{1})
    ishift = cellfun(@(i,d) i+d(p), ind, D, 'unif', 0);
    B(p,dots{:}) = reshape(A(ishift{:})-Amid,[1 m]);
end
B((prod(w)+1)/2,dots{:}) = reshape(Amid,[1 m]);
B = reshape(B, [w m]);

% Bruno

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <63272f32-03e1-455f-aae3-200b8771bec3@googlegroups.com>...
> Well you didn't use blockproc correctly. If you'd like one of my demos using it, let me know.

Please do, it would be great to learn how to use blockproc properly as it seems quite useful. So if I were to call it correctly, it would run in a fraction of the time? I tested it using the Parallel Computing Toolbox, and it dropped down to around 4minutes, but I am guessing I can get those times if I call the function properly...

At Bruno:

Thank you very much for your ingenious solutions. Just wondering, if I know the index of the voxels I want to perform the sliding operation on, is it possible to include this in your solution?

"K M" wrote in message <jum94l$i4k$1@newscl01ah.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <63272f32-03e1-455f-aae3-200b8771bec3@googlegroups.com>...
>
> At Bruno:
>
> Thank you very much for your ingenious solutions. Just wondering, if I know the index of the voxels I want to perform the sliding operation on, is it possible to include this in your solution?

I don't see what why it's not possible.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jumnhv$j5o$1@newscl01ah.mathworks.com>...
> "K M" wrote in message <jum94l$i4k$1@newscl01ah.mathworks.com>...
> > ImageAnalyst <imageanalyst@mailinator.com> wrote in message <63272f32-03e1-455f-aae3-200b8771bec3@googlegroups.com>...
> >
> > At Bruno:
> >
> > Thank you very much for your ingenious solutions. Just wondering, if I know the index of the voxels I want to perform the sliding operation on, is it possible to include this in your solution?
>
> I don't see what why it's not possible.
>
> Bruno

% Data
A=ceil(10*rand(6,8))
w = [3 3];
% each row is the indices of one selected pixel
% it must be inside A with a margin of (w-1)
pixelind = [4 3;
            2 2;
            2 7;
            5 4]

% Engine
n = size(A);
npixels = size(pixelind,1);
dfun = @(w) linspace(1,w,w);
d = arrayfun(dfun, w, 'unif', 0);
D = cell(size(d));
[D{:}] = ndgrid(d{:});
pixellin = num2cell(pixelind,1);
pixellin = sub2ind(size(A),pixellin{:});
D = sub2ind(size(A),D{:});
imid = (prod(w)+1)/2;
B = bsxfun(@plus, D(:)-D(imid), pixellin');
B = A(B);
Amid = B(imid,:);
B = bsxfun(@minus, B, Amid);
B(imid,:) = Amid;
B = reshape(B, [w npixels])

% Bruno