How can I zero out weak elements near stronger ones in a matrix?

Question

I have a pixel matrix containing some points and a lot of zero elements. From those non-zero points, I want to discard those that have a stronger point in range N withing the matrix. The range is an euclidean distance between the pixels.

input  = [0.0 0.0 0.0 0.9 0.0 0.0 
          0.0 0.0 0.2 0.0 0.0 0.5 
          0.0 0.0 0.7 0.0 0.0 0.0 
          0.0 0.4 0.1 0.0 0.0 0.0];

output = [0.0 0.0 0.0 0.9 0.0 0.0   % 0.7 is the largest number in range
          0.0 0.0 0.0 0.0 0.0 0.5   % 0.2 got removed; was next to 0.9 and 0.7
          0.0 0.0 0.7 0.0 0.0 0.0   % 0.7 is the largest number in range
          0.0 0.0 0.0 0.0 0.0 0.0]; % 0.1 and 0.4 both got removed; were next to 0.7

Update: This is what I've came up so far. It iterates over all non-zero pixels and compares the current pixel to the largest of the neighborhood. However, the neighborhood contains to many pixels. Rather than selecting the area through index offsets, I need to select a circular area somehow. Moreover, I'd be grateful if there is a shorter way, maybe replacing the loops with built-in Matlab functions like conv2 or imfilter.

% Discard points near stronger points
points = find(Image > 0);
radius = args.Results.min_dist;
for i = 1:size(points)
    [index_y, index_x] = ind2sub(size(Image), points(i));
    %  Find neighborhood
    from_x = max(index_x-radius, 1);
    from_y = max(index_y-radius, 1);
    to_x = min(index_x+radius, size(Image, 2));
    to_y = min(index_y+radius, size(Image, 1));
    neighbors = Image(from_y:to_y, from_x:to_x);
    % Discard if there is a stronger neighbor
    largest = max(max(neighbors));
    if Image(index_y, index_x) < largest
        Image(index_y, index_x) = 0;
    end
end

Answer 1

If you had the image processing toolbox installed (or one of the free alternatives mentioned in Try Hard's answer), you could make this happen quite easily:

The function imdilate is basically a sliding max filter. So what we do is, for each pixel we look for the maximum in the neighborhood specified by your radius R. Then we compare the actual value with the found maximum and if it is smaller, we set the value to zero.

function A = zeroOutWeakElements(A, R)
[X,Y] = ndgrid(-ceil(R):ceil(R));
neighborhood = (X.^2 + Y.^2)<=R^2;
A(imdilate(A,neighborhood)>A) = 0;

For large full matrices and small distances this would also be significantly faster than the currently accepted solution. The benefit will wear off with sparse matrices and large radii, but I guess you should test with actual data to be sure what's best.

Answer 2

The basic workflow to solve this problem could be described as discussed below -

1. Sort all the non-zero points.
2. Start with the highest valued point and set all non-zeros points within its N-neighbourhood to zeros.
3. Do the same for the second highest one and set all non-zeros points within its N-neighbourhood to zeros excluding the non-zero points with higher value than the current one. This excluding part could be implemented with a very useful MATLAB tool triu inside the codes.
4. Continue till every non-zero point is covered. Of course, as we move down this ladder, there will be lesser points to search for, because of the excluding clause discussed earlier.

These steps could be implemented with a vectorized approach, using no special toolboxes and assuming A as the input -

%// Distance parameter
N = 2; 

%// Find all non-zero points and then sort them in descending manner
[x,y] = find(A~=0)
pts = [x y]
[val,sorted_idx] = sort(A(A~=0),'descend')
pts = pts(sorted_idx,:)

%// Find euclidean distances
distmat = sqrt(squared_dist(pts,pts))

%// Find points to be removed (set to zero); then calculate their linear indices
rm_pts = pts(any(triu(distmat<N,1),1),:)
rm_lin_idx = sub2ind(size(A),rm_pts(:,1),rm_pts(:,2))

%// Use those linear indices to set those in the input as zeros
out = A;
out(rm_lin_idx) = 0;

Associated function code (to find squared euclidean distances) -

function sq_distmat = squared_dist(A,B)

[nA,dim] = size(A);
nB = size(B,1);

A_ext = ones(nA,dim*3);
A_ext(:,2:3:end) = -2*A;
A_ext(:,3:3:end) = A.^2;

B_ext = ones(nB,dim*3);
B_ext(:,1:3:end) = B.^2;
B_ext(:,2:3:end) = B;

sq_distmat = A_ext * B_ext.';

return;

Code run -

A =
         0         0         0    0.9000         0         0
         0         0    0.2000         0         0    0.5000
         0         0    0.7000         0         0         0
         0    0.4000    0.1000         0         0         0
out =
         0         0         0    0.9000         0         0
         0         0         0         0         0    0.5000
         0         0    0.7000         0         0         0
         0         0         0         0         0         0

Answer 3

There are various freely available substitutes for MATLAB image processing toolbox functions; for instance you can look into octave, or if that does not suit you see e.g. this source although it requires mex compilation.

You are looking for a filter with an abrupt cutoff. Filters are a standard tool in the toolbox. Start by looking at the documentation of tools like imfilter and perhaps imdilate and imopen.