How to sum values in a matrix which are in the specified distance from one value?

I have a matrix 1000x1000. It was originally a tiff file; a map in which every pixel was 1km x 1km. It this matrix I have 10 specific values (I know in which columns and rows they are, I also know their indexes in a matrix; those values (previously pixels) are actually localizations of some cities). I need to sum the values from the matrix which are in a distance to those cities no bigger than 50 km. Example: matrix contains the values of temperature; I want to know mean temperature around the cities, 50 km from the middle of each city (it doesn matter if it's sum or mean value; my problem is how to find those 50 km distanced pixels). How can I do this?

Answer 1

[Image Analyst]

What I'd do is to create a mask of circles of radius 50 around the cities. Then use that in a call to regionprops() asking for the mean temperature. Try this:

% Initialization Steps.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 14;

%--------------------------------------------------------------------------------------------
% Read in image -- the temperature map array.
grayImage = imread('concordorthophoto.png');
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels == 3
    grayImage = rgb2gray(numberOfColorChannels);
end
subplot(2, 1, 1);
imshow(grayImage, [])
axis('on', 'image')

%--------------------------------------------------------------------------------------------
% Get random x and y points to indicate city locations, 
% or else read in the ones that you have.
numCities = 10;
x = 50 + (columns - 100) * rand(numCities, 1);
y = 50 + (rows - 100) * rand(numCities, 1);
hold on;
plot(x, y, 'r.', 'MarkerSize', 40)
hold off;
caption = sprintf('Original Image with %d cities marked.', numCities);
title(caption, 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------
% Get a mask with those cities and a radius of 50
% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
[columnsInImage, rowsInImage] = meshgrid(1:columns, 1:rows);
radius = 50;
% Next create the circle in the image.
allCircles = false(rows, columns);
for k = 1 : numCities
    centerX = x(k);
    centerY = y(k);
    circlePixels = (rowsInImage - centerY).^2 ...
        + (columnsInImage - centerX).^2 <= radius.^2;
    % circlePixels is a 2D "logical" array.
    allCircles = allCircles | circlePixels;
end
% Now, display it.
subplot(2, 1, 2);
imshow(allCircles) ;
caption = sprintf('Binary image of %d circles of radius %.1f pixels.', numCities, radius);
title(caption, 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------
% Get the temperatures of those 10 cities and their locations
props = regionprops(allCircles, grayImage, 'MeanIntensity', 'Centroid');
caption = sprintf('Binary image of %d circles of radius %.1f pixels.', length(props), radius);
title(caption, 'FontSize', fontSize);
xy = vertcat(props.Centroid)
x = xy(:, 1);
y = xy(:, 2);
% Get temperatures.  By default, they are sorted from leftmost city to rightmost city.
meanTemps = [props.MeanIntensity]
% Plot temp over city
for k = 1 : length(meanTemps)
    str = sprintf('   City #%d had T = %.1f', k, meanTemps(k));
    text(x(k), y(k), str, 'Color', 'r', 'FontSize', 10, 'FontWeight','bold', 'VerticalAlignment','middle');
end
g = gcf;
g.WindowState = 'maximized';

[Izabela Wojciechowska]

Thank you! Works great. Creating such mask (with cities of specified radius) was exactly what I needed.

Answer 2

[Simon Chan]

You may try the following where the mask indicates the pixel where its distance is 50 pixels from your preset center pixel.

if true
  % clear; clc;
separation = 1;                             % Angle separation, 1 degree
rho = 50;                                   % Radius, 50 pixels
num_dot = 360/separation +1;                % Total number of points
theta = linspace(-pi,pi,num_dot);           % Theta separation
cx = 500;                                    % Your center point, x-coordinates = 500 (example)
cy = 250;                                    % Your center point, y-coordinates = 250 (example)
[xA,yA] = pol2cart(theta,rho);              % Use pol2cart
A = zeros(1000,1000);                         % Original region, example, size: 1000x1000
[Ny, Nx]= size(A);
xA_convert = round(cx-xA);                  % Convert back from polar coordinates
yA_convert = round(cy-yA);
for k = 1:num_dot
  A(yA_convert(k),xA_convert(k)) =1;      % Put the dots in matrix A
end
mask = false(Ny,Nx);
mask(A>0)=true;
figure(1)
plot(cx,cy,'ro');                           % Your center point
hold on
plot(xA_convert,yA_convert,'b+')
xlim([0 Nx])
ylim([0 Ny])
end

[Izabela Wojciechowska]

Thank you for the response.

But where in your code is the sum of values in the matrix inside the circle? It looks like your code only creates the coordinates (columns and rows) of the circle points, am I right?