Plotting random points within Boundary

Hello everyone,

I could use your help to improve the following code. Its about creating a map with random positions of points.

How to put points only inside of country borders?

 C = get(handles.number_cities, 'String');
        C = str2double(C);
        if(isnan(C)), errordlg('Parameters must be numerical.','Mapping Error');
        elseif(C <= 0), errordlg('Number of cities must be greater than zero','Mapping Error');
        else
            M = zeros(C,2);             % Map initialization
            for id_C = 1:1:C
                M(id_C,1) = randi(100); % Random map
                M(id_C,2) = randi(100);
            end

Answer 1

[Image Analyst]

Azra, try this:

% Code to find the boundaries of a country on a map and place a specified number of points in it.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%===============================================================================
% Read in gray scale image.
folder = pwd;
baseFileName = 'map.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
rgbImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(rgbImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
	grayImage = rgbImage(:, :, 1); % Take red channel.
else
	grayImage = rgbImage;
end
% Display the original image.
subplot(2, 2, 1);
imshow(rgbImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;

%===============================================================================
% Segment the image.
% Display the histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Cropped Image Histogram', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);

% Specify a threshold.  The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 240;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);

% Create a binary image
mask = grayImage > threshold;
subplot(2, 2, 3);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;

% Get rid of any holes.
mask = imfill(mask, 'holes');
% Take the largest region.
mask = bwareafilt(mask, 1);
subplot(2, 2, 3);
imshow(mask);
title('Final Mask', 'FontSize', fontSize);
drawnow;

%===============================================================================
% Get the perimeter.
boundaries = bwboundaries(mask);
boundaries = boundaries{1}; % Extract from cell.
xb = boundaries(:, 2);
yb = boundaries(:, 1);

% Get random points.
numPoints = nnz(grayImage); % Get way more points than we need.
% Say how many points we need inside the boundary.
numPointsNeeded = 300;
xp = columns * rand(1, numPoints);
yp = rows * rand(1, numPoints); % You can round them to integers if you want to or need to, like if you're going to use them as indexes in an array.
subplot(2, 2, 4);
imshow(rgbImage);
drawnow;

% Find out which points are in the polygon.
numPointsInside = 0;
insideBoundary = false(numel(xp), 1);
for k = 1 : numel(xp)
	% Check if (xp, yp) is inside the boundaries defined by xb and yb.
	insideBoundary(k) = inpolygon(xp(k), yp(k), xb, yb);
	if insideBoundary(k)
		numPointsInside = numPointsInside + 1;
		if numPointsInside >= numPointsNeeded
			break;
		end
	end
end
% Plot them all
x = xp(insideBoundary);
y = yp(insideBoundary);
hold on;
plot(x, y, 'r.', 'MarkerSize', 10);
caption = sprintf('Original image with %d locations', numPointsInside);
title(caption, 'FontSize', fontSize);

msgbox('Done!  Thank you Image Analyst!');

[Image Analyst]

I edited the above to correct a small mistake where I had rows and columns swapped whn I created xp and yp. Anyway, here is another demo that is similar where I used the RGB-to-HSV transform to find the country and then also found the centroid and plotted it as a blue crosshairs over the image.

% Code to find the boundaries of a country on a map and place a specified number of points in it.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%===============================================================================
% Read in gray scale image.
folder = pwd;
baseFileName = 'random_points_USA_Map.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
rgbImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(rgbImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
	grayImage = rgbImage(:, :, 1); % Take red channel.
else
	grayImage = rgbImage;
end
% Display the original image.
subplot(2, 2, 1);
imshow(rgbImage);
axis('on', 'image');
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : "%s"', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;

%===============================================================================
% Segment the image.
% Convert to HSV Color Space.
hsvImage = rgb2hsv(rgbImage);
grayImage = hsvImage(:,:,2); % Take saturation channel to find vividly colored regions.
h1 = subplot(2, 2, 2);
imshow(grayImage, []);
axis('on', 'image');
impixelinfo
grid on;
title('Saturation Image', 'FontSize', fontSize);

% Specify a threshold for the saturation image.  The image is double in the range 0-1 so the threshold will be a fraction.
threshold = 0.05;

% Create a binary image
mask = grayImage > threshold;
h2 = subplot(2, 2, 3);
imshow(mask);
title('Initial Mask', 'FontSize', fontSize);
drawnow;

% Get rid of any holes.
mask = imfill(mask, 'holes');
% Take the largest region.
mask = bwareafilt(mask, 1);
h3 = subplot(2, 2, 3);
imshow(mask);
axis('on', 'image');
hold on;
title('Final Mask.  Centroid at Blue +', 'FontSize', fontSize);
drawnow;

%==============================================================================================================================================
% Get the perimeter.
boundaries = bwboundaries(mask);
boundaries = boundaries{1}; % Extract from cell.
xb = boundaries(:, 2);
yb = boundaries(:, 1);

% Get random points.
numPoints = nnz(grayImage); % Get way more points than we need.
% Say how many points we need inside the boundary.
numPointsNeeded = 800;
xp = columns * rand(1, numPoints);
yp = rows * rand(1, numPoints); % You can round them to integers if you want to or need to, like if you're going to use them as indexes in an array.
h4 = subplot(2, 2, 4);
imshow(rgbImage);
axis('on', 'image');
hold on;
drawnow;

% Find out which points are in the polygon.
numPointsInside = 0;
insideBoundary = false(numel(xp), 1);
for k = 1 : numel(xp)
	% Check if (xp, yp) is inside the boundaries defined by xb and yb.
	insideBoundary(k) = inpolygon(xp(k), yp(k), xb, yb);
	if insideBoundary(k)
		numPointsInside = numPointsInside + 1;
		if numPointsInside >= numPointsNeeded
			break;
		end
	end
end
% Plot them all
x = xp(insideBoundary);
y = yp(insideBoundary);
plot(h3, x, y, 'r.', 'MarkerSize', 10);
plot(h4, x, y, 'r.', 'MarkerSize', 10);
caption = sprintf('Original image with %d locations', numPointsInside);
title(caption, 'FontSize', fontSize);

% Just for fun, let's find the geographic centroid.
props = regionprops(mask, 'Centroid');
xCentroid = props.Centroid(1);
yCentroid = props.Centroid(2);
plot(h3, xCentroid, yCentroid, 'b+', 'MarkerSize', 50, 'LineWidth', 2);
plot(h4, xCentroid, yCentroid, 'b+', 'MarkerSize', 50, 'LineWidth', 2);

msgbox('Done!  Thank you Image Analyst!');

[Adam Danz]

+1 for cool demos

[azra hosic]

When I run the code, I get this:

No public property WindowState exists for class matlab.ui.Figure.

Error in random_points_USA_map (line 46)

g.WindowState = 'maximized';

Answer 2

[Adam Danz]

If you have the perimeter coordinates, you can create a bunch of random numbers and then eliminate the ones outside of the perimeter using inpolygon.

If you don't have the perimeter coordinates, one way of getting them is with the borders function on the file exchange.

Example:

clf()
[lat,lon] = borders('Bosnia and Herzegovina');
plot(lon,lat,'k-')
axis equal
hold on
grid on

% Produce n random coordinates within the
% extent of latitude and longitude
n = 100;
xy = inf(n,2);
in = false;
while ~all(in)
    % Determine which random coordinates are within
    % the perimeter
    in = inpolygon(xy(:,1),xy(:,2),lon,lat);
    
    % Replace external coordinates with new rand values
    xy(~in,1) = rand(sum(~in),1) * range(lon) + min(lon);
    xy(~in,2) = rand(sum(~in),1) * range(lat) + min(lat);
end

% Add random coordinates
plot(xy(:,1), xy(:,2), 'rx')

xlabel('lon')
ylabel('lat')
title(char([1041 1086 1089 1085 1072 32 1080 32 1061 1077 1088 1094 1077 1075 1086 1074 1080 1085 1072]))

[azra hosic]

Also here, when I run the code I get this message:

Undefined function or variable 'borders'.

Error in karta (line 2)

[lat,lon] = borders('Bosnia and Herzegovina');

[Adam Danz]

Check out the link in the second sentence of my answer.