Reputation: 117
This is the following part of the below:
2) Additional question:
After getting the average of the non-zero neighbors, I also want to test if the neighbor elements are equal, lesser, or greater than the average of the nonzeros. If it is greater or equal then '1' or else '0'.
Note: if the neighbors are with in the radius of the two or more centres, take the smallest centre average to test.
0 12 9
4 **9** 15
11 19 0
The '9' in the middle is within the radius of 12, 15, and 19 centres, so take the minimum average of those min[9.000, 9.000, 8.000]=8.000
For example, when radius = 1 m or 1 element away.
new_x =
0 0 0 0 0
0 0 **9.0000** 9.0000 0
0 4.0000 9.0000 **9.0000** 0
0 **8.3333** **8.0000** 0 0
0 2.0000 4.0000 8.0000 0
0 4.0000 5.0000 8.0000 0
0 0 0 0 0
Test_x =
0 0 0 0 0
0 0 **9.0000** 1 0
0 0 1 **9.0000** 0
0 **8.3333** **8.0000** 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
=================================================================================
1) Say if I have a matrix, shown as below,
X =
0 0 0 0 0
0 0 12 9 0
0 4 9 15 0
0 11 19 0 0
0 2 4 8 0
0 4 5 8 0
0 0 0 0 0
and I want to find the average of the surrounding non-zero elements that is greater than 10. The rest of the elements still remain the same i.e. elements < 10.
So I want my solution to look something like,
new_x =
0 0 0 0 0
0 0 9.0000 9.0000 0
0 4.0000 9.0000 9.0000 0
0 8.3333 8.0000 0 0
0 2.0000 4.0000 8.0000 0
0 4.0000 5.0000 8.0000 0
0 0 0 0 0
Not: that I am NOT only looking at the neighbors of the element thats greather than some value (i.e. 10 in this case).
Lets say any elements thats are greater than 10 are the 'centre' and we want to find the avearge of the non-zeros with the radius of say 1 m. where 1 metre = 1 element away from the centre.
Note: It might not always be 1 meter away in radius i.e. can be 2 or more. In this case it wont be just top, bottom, left and right of the centre.
****Also Be aware of the matrix boundary. For example, when radius = 2 or more, some of the average of nonzero neighbors are out side the boundary.**
For example,
For radius =1 m = 1 element away, new_x = average of [(i+1,j) , (i-1,j) , (i,j+1) and (i,j-1)] - top, bottom, right, and left of the centre.
For radius =2 m = 2 elements away, new_x = average of [(i+1,j), (i+2,j) , (i-1,j) , (i-2,j), (i,j+1), (i,j+2), (i,j-1), (i,j-2), (i+1,j+1), (i+1,j-1), (i-1,j-1), and (i-1,j+1)].
==================================================================
I have tried a few things before, however I am not familiar with the functions.
So please help me to solve the problem.
Thank you in advance.
Upvotes: 5
Views: 2439
Reputation: 124563
EDIT: Note this requires functions from the Image Processing Toolbox, namely: COLFILT and STREL
r = 1; %# radius
t = 10; %# threshold value
mid = round((2*r+1)^2/2); %# mid point
nhood = getnhood(strel('diamond', r));
nhood(mid) = false;
fcn = @(M)sum(M(nhood(:),:),1)./(sum(M(nhood(:),:)~=0)+all(M(nhood(:),:)==0)).*(M(mid,:)>=t)+M(mid,:).*(M(mid,:)<t);
new_x = colfilt(x, 2*[r r]+1, 'sliding',fcn)
For r=1:
new_x =
0 0 0 0 0
0 0 9 9 0
0 4 9 9 0
0 8.3333 8 0 0
0 2 4 8 0
0 4 5 8 0
0 0 0 0 0
For r=2:
new_x =
0 0 0 0 0
0 0 11.2 9 0
0 4 9 10.167 0
0 7 7.7778 0 0
0 2 4 8 0
0 4 5 8 0
0 0 0 0 0
In fact, it should work for any radius >= 1
Notice how the diamond shape structuring element represents the neighborhood:
nhood =
0 1 0
1 0 1
0 1 0
nhood =
0 0 1 0 0
0 1 1 1 0
1 1 0 1 1
0 1 1 1 0
0 0 1 0 0
and so on..
Explanation:
We use the COLFILT function which traverse the matrix using a sliding neighborhood of NxN, and places each block as a column in a temporary matrix.
We process each column of this temp matrix (blocks) using the function fcn
, and the result will be placed in the correct location once finished (COLFILT uses IM2COL and COL2IM underneath).
We check for two cases depending of the value of the center of the block:
If its less than 10, it returns that value unchanged: M(mid,:)
if its >=10, we compute the mean of the non-zero elements of its neighborhood
sum(M(nhood(:),:),1) ./ (sum(M(nhood(:),:)~=0) + all(M(nhood(:),:)==0))
.
The last term in there is necessary to avoid dividing by zero
Notice how the result of 1 & 2 above are combined using R1.*(M(mid,:)<t) + R2.*(M(mid,:)>=t)
to emulate an if/else choice.
Upvotes: 4
Reputation: 1193
EDIT: Here is a solution that does not require the Image Processing Toolbox. It does, however, use conv2nan.m
which is part of the free NaN toolbox.
This approach relies on doing two different filtering/convolution operations: one that gets the sum of surrounders for each element, and one that gets the count of nonzero surrounders. Then, you are ready to combine them to get the average of nonzero surrounders only. Like this:
% set up starting point matrix with some zeros
X = magic(4);
X(X < 5) = 0;
X(X == 0) = NaN; % convert zeros to NaNs to work with conv2nan
countmat = double(X > 0);
cmat = [0 1 0;
1 0 1;
0 1 0]; % consider surrounding elements only
[m1,c] = conv2nan(X,cmat); % sum of surrounding elements
[m2,c] = conv2nan(countmat,cmat); % number of surrounding elements > 0
x_new = m1./m2; % compute average we want
x_new = x_new(2:end-1,2:end-1); % trim edges created by conv2
x_new(~countmat) = 0; % restore zero elements
x_new(X < 10) = X(X < 10) % restore < 10 elements
It does some extra work in that the convolutions are done for all elements and not just those that are >= 10. But it's more general than the manual looping approach.
Upvotes: 1
Reputation: 125854
Here is the algorithm I think you are describing in your question. For each pixel:
If this is correct (as it appears to be from the sample matrices you gave), then you could use the function NLFILTER from the Image Processing Toolbox (if you have access to it) to perform this operation:
fcn = @(x) [x(5) sum(x(2:2:8))/max(sum(x(2:2:8) > 0),1)]*[x(5) < 10; x(5) >= 10];
new_x = nlfilter(X,[3 3],fcn);
EDIT: If you don't have access to the Image Processing Toolbox, you can also do this using the built-in CONV2 function, like so:
kernel = [0 1 0; ... %# Convolution kernel
1 0 1; ...
0 1 0];
sumX = conv2(X,kernel,'same'); %# Compute the sum of neighbors
%# for each pixel
nX = conv2(double(X > 0),kernel,'same'); %# Compute the number of non-zero
%# neighbors for each pixel
index = (X >= 10); %# Find logical index of pixels >= 10
new_x = X; %# Initialize new_x
new_x(index) = sumX(index)./max(nX(index),1); %# Replace the pixels in index
%# with the average of their
%# non-zero neighbors
The above handles your radius = 1 case. To address your radius = 2 case, you just have to change the convolution kernel to the following and rerun the above code:
kernel = [0 0 1 0 0; ...
0 1 1 1 0; ...
1 1 0 1 1; ...
0 1 1 1 0; ...
0 0 1 0 0];
Upvotes: 3
Reputation: 189676
You could do something like this: (tested in Octave, should work in matlab)
octave-3.2.3:17> toohigh = (x>=10)
toohigh =
0 0 0 0 0
0 0 1 0 0
0 0 0 1 0
0 1 1 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
octave-3.2.3:18> nbr_avg = filter2(ones(3,3)/9,x)
nbr_avg =
0.00000 1.33333 2.33333 2.33333 1.00000
0.44444 2.77778 5.44444 5.00000 2.66667
1.66667 6.11111 8.77778 7.11111 2.66667
1.88889 5.44444 8.00000 6.11111 2.55556
1.88889 5.00000 6.77778 4.88889 1.77778
0.66667 1.66667 3.44444 2.77778 1.77778
0.44444 1.00000 1.88889 1.44444 0.88889
octave-3.2.3:19> y=x; y(toohigh) = nbr_avg(toohigh)
y =
0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 5.44444 9.00000 0.00000
0.00000 4.00000 9.00000 7.11111 0.00000
0.00000 5.44444 8.00000 0.00000 0.00000
0.00000 2.00000 4.00000 8.00000 0.00000
0.00000 4.00000 5.00000 8.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000
The filter2
function allows you to filter on neighbors (not sure what function you want...), and if you use a boolean index matrix (toohigh
in this case) to select those members of the original matrix that are too high, you can replace them with the ones you want.
More specifically, filter2 allows you to convolve with an arbitrary matrix. The matrix of all ones does a spatial low pass filter.
note: my math doesn't match yours. I'm not quite sure why you want to average only the nonzero neighbors (that gives higher weight to nonzero neighbors when there are zeros), but if you wanted to do that, you could do filter2(ones(3,3),x) ./ M
where M = filter2(ones(3,3),(x ~= 0))
is the count of nonzero neighbors.
Upvotes: 1