# Thread Subject: Difference between pixel values and centre pixel in a local neighbourhood

 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 20 Jul, 2012 10:55:17 Message: 1 of 13 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 !
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: ImageAnalyst Date: 21 Jul, 2012 01:40:13 Message: 2 of 13 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);
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 21 Jul, 2012 18:08:16 Message: 3 of 13 ImageAnalyst 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
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: ImageAnalyst Date: 21 Jul, 2012 22:48:21 Message: 4 of 13 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.
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 23 Jul, 2012 16:20:24 Message: 5 of 13 ImageAnalyst wrote in message ... > 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?
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 23 Jul, 2012 17:08:09 Message: 6 of 13 "K M" wrote in message ... > ImageAnalyst wrote in message ... > > 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)
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 23 Jul, 2012 20:20:10 Message: 7 of 13 ImageAnalyst wrote in message ... > 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.
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: Bruno Luong Date: 23 Jul, 2012 21:52:09 Message: 8 of 13 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
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: ImageAnalyst Date: 24 Jul, 2012 00:38:21 Message: 9 of 13 Well you didn't use blockproc correctly. If you'd like one of my demos using it, let me know.
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: Bruno Luong Date: 24 Jul, 2012 08:01:13 Message: 10 of 13 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
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: K M Date: 24 Jul, 2012 13:49:09 Message: 11 of 13 ImageAnalyst 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?
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: Bruno Luong Date: 24 Jul, 2012 17:55:11 Message: 12 of 13 "K M" wrote in message ... > ImageAnalyst 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
 Subject: Difference between pixel values and centre pixel in a local neighbourhood From: Bruno Luong Date: 25 Jul, 2012 06:58:09 Message: 13 of 13 "Bruno Luong" wrote in message ... > "K M" wrote in message ... > > ImageAnalyst 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

### Everyone's Tags:

Separated by commas
Ex.: root locus, bode

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.