Take 3x3 nearest neighbours round centroids of all spots and find the median of their mean

Hello, I have an image (TIFF) that I calculate the centroids and extract the median intensity of all these centroids.

I have the centroid positons in x1,y1 and loop through them the obtin the intensity of each spot and then take the median

Intensity(i)=IM(y1(i),x1(i));   %Centroid intensity
...
medIntensity=median(Intensity(:));

I now want to take a 3x3 nearest neighbours round each of the x1,y1 points and obtain the median of their mean.

Are there any suggestions how I achieve this.

Thanks

Jason

Answer 1

[Matt J]

sz=size(IM);

centroidIdx=sub2ind(sz,y1,x1);

T=combinations(1:3,1:3); 
delta=sub2ind(sz,T{:,2},T{:,1})'-(sz(1)+2);

result=median(  mean( IM(centroidIdx(:) + delta),2) )

[Jason]

Thanks Matt for both your answers. I hadn't heard of the combinations function before, but I prefer this as its simpler.

How would I go about plotting the locations of all these coordiantes (i.e. the 3x3 square) actually on my image?

This is my attempt that doesnt work

      sz=size(IM);
        
      centroidIdx=sub2ind(sz,round(y1r),round(x1r))   %x1r, y1r and not integers so round
      T=combinations(1:3,1:3); 
      delta=sub2ind(sz,T{:,2},T{:,1})'-(sz(1)+2);
      result=median(  mean( IM(centroidIdx(:) + delta),2) )
      ReportMessage(app,['Median(3x3 average) = ',num2str(result)])   
      A=centroidIdx(:) + delta

      [row,col] = ind2sub(IM,A)
      hold(ax,'on');
      plot(ax,col,row,'g.')

[Matt J]

  IM=rand(12,14);
  C=combinations(2:4:12,2:5:14);[y1r,x1r]=deal(C{:,1}, C{:,2});


      sz=size(IM);
        
      centroidIdx=sub2ind(sz,round(y1r),round(x1r));   %x1r, y1r and not integers so round
      T=combinations(1:3,1:3); 
      delta=sub2ind(sz,T{:,2},T{:,1})'-(sz(1)+2);

      A=centroidIdx(:) + delta;

      [Y,X] = ind2sub(sz,A);
      
      imagesc(IM);colormap(gray); hold on
      plot(x1r,y1r,'m+',X(:),Y(:),'co','MarkerSize',11); hold off

●

[Jason]

Perfect, thanks

[Jason]

Hi matt, Im revisiting this - how woudl this be modified to allow e.g. a 5x5 or even 6x6 or generally nxn instead of just 3x3?

[Matt J]

  IM=rand(24,28);
  n=5;  
  hn=ceil(n/2);
  
  C=combinations(hn:6:height(IM),hn:8:width(IM));  
  [y1r,x1r]=deal(C{:,1}, C{:,2});


      sz=size(IM);
        
      centroidIdx=sub2ind(sz,round(y1r),round(x1r));   %x1r, y1r and not integers so round
      
      
      T=combinations(1:n,1:n); 
      delta=sub2ind(sz,T{:,2},T{:,1})'- sub2ind(sz,hn,hn);

      A=centroidIdx(:) + delta;

      [Y,X] = ind2sub(sz,A);
      
      imagesc(IM);colormap(gray); hold on
      plot(x1r,y1r,'m+',X(:),Y(:),'co','MarkerSize',11); hold off

[Jason]

Hi thanks for this, but now this part doesn't work

 result=median( mean(IM(centroidIdx(:) + delta),2) )

[Matt J]

It should work, but be mindfult that your centroids must be deep enough in the interior of the image that the whole of each nxn neighborhood fits within the image boundaries.

  IM=rand(24,30);
  n=5;  
  hn=ceil(n/2);
  
  C=combinations(hn:6:height(IM),hn:8:width(IM));  
  [y1r,x1r]=deal(C{:,1}, C{:,2});


      sz=size(IM);
        
      centroidIdx=sub2ind(sz,round(y1r),round(x1r));   %x1r, y1r and not integers so round
      
      
      T=combinations(1:n,1:n); 
      delta=sub2ind(sz,T{:,2},T{:,1})'- sub2ind(sz,hn,hn);

   result=median( mean(IM(centroidIdx(:) + delta),2) )

result = 0.5203

Answer 2

[Matt J]

movingMean=conv2(IM,ones(3)/9,'same');

centroidIdx=sub2ind(size(IM),y1,x1);

result=median(movingMean(centroidIdx))

Answer 3

[Image Analyst]

If you've used regionprops to find the centroids or the spots, then why not simply also ask it for the mean intensities?

props = regionprops(mask, 'Centroid', 'MeanIntensity');
xyCentroids = vertcat(props.Centroid);
meanIntensities = [props.MeanIntensity]
medianMeanIntensity = median(meanIntensities)

If you want different diameters around the spots, simply call imdilate to enlarge the spots mask.

[Jason]

Thanks IA, Im actually using normalised cross correlation to get the centroids normxcorr2(template,IM);

[Image Analyst]

If you already know in advance the (row, column) or (x,y) coordinates where you want to look, then you can just use conv2 to get the mean in windows around them. For example,

windowSize = 3; % Adjust as needed to other sizes.
kernel = ones(windowSize)/windowSize^2; % You can modify this for larger circles if you don't want square shapes.
sumImage = conv2(grayImage, kernel, 'same');

If you know that, say, the location you want to look at is row 40, column 50 (in other words (x,y) = (50, 40) because x,y is switched from row,column) then you can do

meanIntensity = sumImage(40, 50);

Or to do in a loop for all locations if you have (x,y) in an array with x (column) in column 1 and y (row) in column 2):

numPoints = size(xyCentroids, 1);
means = zeros(numPoints, 1);
for k = 1 : numPoints
    row = xyCentroids(k, 2); % y
    col = xyCentroids(k, 1); % x
    means(k) = sumImage(row, col);
    fprintf('The mean at row=%d, col=%d is %f.\n', row, col, means(k))
end
medianOfMeans = median(means)

[Jason]

Thanks IA, this is awesome. I like to visually see what its actually making measurements on. Is there anyway to get the x,y locations of the kernel so I can plot these onto the image also - again just to confirm visually?

[Image Analyst]

numPoints = size(xyCentroids, 1);
means = zeros(numPoints, 1);
hold on;
for k = 1 : numPoints
    row = xyCentroids(k, 2); % y
    col = xyCentroids(k, 1); % x
    % Overlay a red square on the centroid.
    plot(col, row, 'rs', 'MarkerSize', 30); % Adjust marker size as needed.
    means(k) = sumImage(row, col);
    fprintf('The mean at row=%d, col=%d is %f.\n', row, col, means(k))
end
medianOfMeans = median(means)

[Jason]

Thanks again IA. Its not really what Im after. I want to extend as below but allow for other sizes besides 3x3. Ill provbabaly need 5x5 or 7x7 only.

I have the centroids (red circles). I then want to take e.g. a 3x3 set of pixels surrounding that centroid and calc the median intensity. I have this using Matt J's answer.

The bit I can't work out is to increase say to 5x5 (and have them plotted as in the image below

    sz=size(IM);        
    centroidIdx=sub2ind(sz,round(y1),round(x1));   %x1, y1 and not integers so round
    T=combinations(1:3,1:3);  % I thought I could change this to (1:5,1:5) would work
    delta=sub2ind(sz,T{:,2},T{:,1})'-(sz(1)+2);
    result=median( mean(IM(centroidIdx(:) + delta),2));
    ReportMessage(app,['Median(3x3 average) = ',num2str(result)]);
    A=centroidIdx(:) + delta;
    [Y,X] = ind2sub(sz,A);      
    hold(ax,'on');
    plot(ax,round(x1),round(y1),'r.',X(:),Y(:),'c.','MarkerSize',11); hold(ax,'off')

So when I try a 5x5, I do this:

    T=combinations(1:5,1:5);    % T=combinations(1:3,1:3); 
    delta=sub2ind(sz,T{:,2},T{:,1})'-(sz(1)+2); 

But it gives this (i.e the 5x5 are not centred to the centroid)

Thanks

Jason

[Image Analyst]

Not sure you understand what conv2 and imfilter do. They scan the image with a window of specified size, and give a local mean at the center of the window location. As I understand it, this is exactly what you want. All you need to do is to ask for the intensity of the output image at the desired locations. And it's far, far less confusing than the combinations() and sub2ind() stuff of the code you're using. There will be no offset away from the center, like you're getting, as long as you use odd-numbered window widths. The only offset will be a faction of a pixel but you seem to be OK with that since you know that the centroids do not fall on integer pixel locations and you need to round the coordinates of the centroid to the nearest pixel location. My code only put a marker at the centroid but I assure you it is taking into account all pixels perfectly centered around the the rounded centroid pixel location. If you want, you could put a rectangle or grid of markers in the overlay around each centroid you probe, but of course that doesn't affect the results.

If you want more accurate means (like to handle centroids off the original pixels) then you'd have to resize the image to a larger image with imresize().

[Jason]

OK, thanks

"If you want, you could put a rectangle or grid of markers in the overlay around each centroid you probe"

Yes this is exactly what I want. so how to do this for the case of a 3x3 or 5x5 grid?

Interesitng point about imresize to get sub pixel resolution, I assume it uses some kind of interpolation?

[Image Analyst]

For a rectangle, you'd put a rectangle up with the rectangle function. Set the position input with [xLeft, yTop, width, height]. So if (x,y) is the coordinate (meaning column, row), in the loop over all your centroids, for a 5x5 rectangle you'd have

rectangle('Position', [x-2.5, y-2.5, 5, 5], 'Color', 'c', 'LineWidth', 1);

but you'd have to see where the centroid is after you round it. Is it in the "middle" of a pixel, or is it in the upper left corner of the pixel? You may have to adjust the xLeft and yTop to make sure the rectangle encompasses the 25 pixels you want to encompass.

If you want small markers, then you could do a loop over the pixels surrounding the current centroid, so in the loop you'd havve a double loop for the markers. Something like this (untested)

windowWidth = 5;% Adjust as needed to other sizes.
halfWindowWidth = floor(windowWidth/2);
% Set up a window (kernel) to get the local means at every pixel location.
kernel = ones(windowWidth)/windowWidth^2; % You can modify this for larger circles if you don't want square shapes.
% Get the local means image.  The value at each pixel will be the
% mean of pixels within the window centered at each pixel location.
% Even though conv2 processes the whole image, it's very fast
% because it's very highly optimized.
sumImage = conv2(grayImage, kernel, 'same');
% Now get the mean within the window for the few number of discrete
% locations that we've gotten prior in xyCentroids.
numPoints = size(xyCentroids, 1);
means = zeros(numPoints, 1);
hold on;
for k = 1 : numPoints
    % Assume xyCentroids have already been rounded to integers.
    % Convert the x,y into row,column 
    % so we can index into the array properly.
    row = xyCentroids(k, 2); % y
    col = xyCentroids(k, 1); % x
    
    % Overlay a red square on the centroid.
    plot(col, row, 'rs', 'MarkerSize', 30); % Adjust marker size as needed.
    
    % Extract the mean within the window for this particular location.
    means(k) = sumImage(row, col);
    fprintf('The mean at row=%d, col=%d is %f.\n', row, col, means(k))
    
    % Place markers over the nxn pixels in the window centered around the
    % rounded centroid location.
    for xm = col - halfWindowWidth : col + halfWindowWidth
        for ym = row - halfWindowWidth : y + halfWindowWidth
            plot(xm, ym, 'c.', 'MarkerSize', 9);
        end
    end
end
% Get the median of all the mean intensities.
medianOfMeans = median(means)

Be careful with normalized cross correlation, because the centroids can be shifted from where you think they should be if you don't do the indexing correctly. It's kind of tricky so I invite you to see how I handled it in my attached demo. Study that and make sure you understand how I find the location of the matching window.

[Jason]

Thankyou, this is exactly what I want

[Jason]

So here is the version as per your suggestion (thankyou so much)

1: Whats the best way to handle where the 5x5 square extends out of the image (i.e. where the green arrow is). I see there are padding options in conv2, but ideally I'd just like to ignore these ones. Perhaps I just go through all my centroids and remove any within 5 pixels of the edge of the image?

2: Rather than means, I'd prefer the median value within each 5x5 - is there a Kernel for this?

3: Also how woudl the Kernel change for a circular shape rather than square