Connecting points to polygon
Show older comments
Hi,
i'm currently struggling with following issue: I have a table geometry_points_outline which stores the outline of a geometry as point coordinates x and y. the geometry consists of three seperate polygons. I'm able to connect the points with the following code section:
for i = 1:height(geometry_points_outline)
nearestPoints = nearestIndices(i, :);
plot([geometry_points_outline.xCoordinates(i), geometry_points_outline.xCoordinates(nearestPoints(2))], ...
[geometry_points_outline.yCoordinates(i), geometry_points_outline.yCoordinates(nearestPoints(2))], 'k-', 'LineWidth', 2);
plot([geometry_points_outline.xCoordinates(i), geometry_points_outline.xCoordinates(nearestPoints(3))], ...
[geometry_points_outline.yCoordinates(i), geometry_points_outline.yCoordinates(nearestPoints(3))], 'k-', 'LineWidth', 2);
end
I have another table which stores regular space grid points. Now i want to do two things:
- filling the space between the lines so the geometry is represented clearly.
- Find the points of the grid which lie between these lines and and delete them from the grid table.
I tried various things, including drawing a polygon manually or selecting the points manually, but the have not been successful. Do you have any idea how to solve this issue?
I attached both tables as csv files.
Thanks for your help!
Accepted Answer
Your code can't be run because nearestIndices is not provided.
Regardless, I would recommend using polyshape to represent and plot the polygons. Likewise isinterior can be used to assess whether points lie inside the polygons.
20 Comments
Matt J
on 13 Jan 2024
Still doesn't run. See above.
Providing the whole code is overkill, however. Just attach a .mat file with what is needed to run your loop.
Debora Baumann
on 14 Jan 2024
Edited: Walter Roberson
on 14 Jan 2024
Here's one way to obtain the polyshapes.
%% Points outline
% Assuming you have a CSV file with coordinate points, load the data
coordinateData = readtable('outline_points_straight.csv'); % Replace with your actual file name
% Extract x and y coordinates from the table
xCoordinates = coordinateData.x;
yCoordinates = coordinateData.y;
L1 = [ 0.056371763743440 0.081642635139480 0.210448634448849]';
L2 = [ 0.009821993289040 0.007773140937910 0.007627949415919]';
XY1=[xCoordinates, yCoordinates, xCoordinates.^0];
G1=XY1*L1<0; %Group points into 3 regions.
G3=XY1*L2>0;
G2=~(G1|G3);
fn=@(G) mkpshape(xCoordinates, yCoordinates,G);
p(1)=fn(G1); %form polyshapes
p(2)=fn(G2);
p(3)=fn(G3);
%Visualize
h=plot(p); hold on
[h.FaceColor]=deal('g','b','m');
[h.EdgeColor]=deal('none');
plot(xCoordinates,yCoordinates,'.'); hold off
axis equal;
axis([1.2280 1.9800 -3.6541 -3.0610])
function p=mkpshape(x,y,G)
%Trace the boundary of the given points and convert to polyshapes
x=x(G); y=y(G);
V=[x,y];
V(end,:)=[];
W=nan(size(V));
v0=V(1,:); V(1,:)=inf;
[d,i]=min(vecnorm(V-v0,2,2));
v1=V(i,:); V(i,:)=inf;
W(1:2,:)=[v0;v1];
for k=3:height(W)
[d,i]=min(vecnorm(V-v1,2,2));
v1=V(i,:); V(i,:)=inf;
W(k,:)=v1;
end
p=polyshape(W,'Simplify',0);
end
Debora Baumann
on 15 Jan 2024
That's amazing! Thank you very much :) How did you determine L1 and L2? I have a very similar geometry i need to visualize too (see attached picture) and I'm guessing that mistake in the plot is do to wrong L1 and L2 values. 
And i need the polygons to be black and i already switched this line [h.FaceColor]=deal('k','k','k'); , but they keep appearing grey.
Thanks!
Matt J
on 15 Jan 2024
And i need the polygons to be black and i already switched this line [h.FaceColor]=deal('k','k','k');, but they keep appearing grey.
Add the line,
[h.FaceAlpha]=deal(1)
How did you determine L1 and L2?
I drew separating lines L1 and L2 on the plot as below using imline. Using the equations for these lines, you can tell whether a point is above or below each of these lines.
Debora Baumann
on 16 Jan 2024
Edited: Debora Baumann
on 16 Jan 2024
I can only draw one line and the endpoints are stored the second i finish my first attempt.
You can do,
h1=imline; wait(h1)
This will cause Matlab to wait until you have finalized the position of the line by double clicking on it.For additional lines, you must make additional calls to imline.
I don't really understand how you have 3 numbers in L1 and L2.
A general line follows the equation A*x+B*y+C=0, which can also be represented L*[x;y;1]=0 where L is the coefficient vector L=[A,B,C]. In this form, you can test whether x,y lies on one side of the line L*[x;y;1]>0 or the other L*[x;y;1]<0.
Note: you can obtain L given any two points on the line (x1,y1) (x2,y2) by doing,
L=cross( [x1,y1,1] , [x2,y2,1])
but keep in mind that L is only unique up to a multiplicative factor. In other words, if L=[A,B,C] is a vector of coefficients for a line, then L=[s*A,s*B,s*C] represents the same line for any non-zero scalar s.
Debora Baumann
on 16 Jan 2024
Debora Baumann
on 17 Jan 2024
Okay, i'm getting there :) I appreciate your help a lot!
Now i don't understand why this is happening:
I plotted my Ls as Lines and it should work and i believe the Gs are written correctly too. Do you have any idea why this is happening?
clear all
%% Points outline
% Assuming you have a CSV file with coordinate points, load the data
coordinateData = readtable('outline_points_straight.csv'); % Replace with your actual file name
xCoordinates = coordinateData.x;
yCoordinates = coordinateData.y;
xCoordinates = coordinateData.x;
coordinateData.x = xCoordinates + 0.363;
xCoordinates = coordinateData.x;
yCoordinates = coordinateData.y;
coordinateData.y = yCoordinates + 4.342;
yCoordinates = coordinateData.y;
%%
endpoints_oZ_straight_left= [1.31455642023346 1.10814785992218;1.99159922178988 0.63927626459144];
endpoints_oZ_straight_right= [1.87795187943582 1.46243518189202;2.23533829024068 0.975363085155054];
%% left line
P1l = [endpoints_oZ_straight_left(1, 1), endpoints_oZ_straight_left(2, 1),1];
%% % Starting and endpoint from line
x_left = [endpoints_oZ_straight_left(1, 1), endpoints_oZ_straight_left(2, 1)];
y_left = [endpoints_oZ_straight_left(1, 2), endpoints_oZ_straight_left(2, 2)];
coefficients_left = polyfit(x_left, y_left, 1);
slope_left = coefficients_left(1);
intercept_left = coefficients_left(2);
% Display the equation of the line
disp(['Equation of the left line: y = ', num2str(slope_left), ' * x + ', num2str(intercept_left)]);
%% left line
% Create a vector with values between 0 and 5 with steps of 0.1
y_values = 0:0.1:5;
% Iterate over each y value
for y_left = y_values
% Solve for x
x_left = (y_left - intercept_left) / (slope_left);
% Calculate Cx and Cy
P2l=[x_left,y_left,1];
L1= cross(P1l,P2l);
% Check if L1(1,1) and L1(1,2)
if L1(1,1) > 0 && L1(1,2) > 0
% If positive, break out of the loop
break;
end
end
%% right line
P1r = [endpoints_oZ_straight_right(1, 1), endpoints_oZ_straight_right(2, 1),1];
%% % Starting and endpoint of line
x_right = [endpoints_oZ_straight_right(1, 1), endpoints_oZ_straight_right(2, 1)];
y_right = [endpoints_oZ_straight_right(1, 2), endpoints_oZ_straight_right(2, 2)];
coefficients_right = polyfit(x_right, y_right, 1);
slope_right = coefficients_right(1);
intercept_right = coefficients_right(2);
disp(['Equation of the right line: y = ', num2str(slope_right), ' * x + ', num2str(intercept_right)]);
%%
% Iterate over each y value
for y_right = y_values
% Solve for x
x_right = (y_right - intercept_right) / (slope_right);
% Calculate Cx and Cy
P2r=[x_right,y_right,1];
L2= cross(P1r,P2r);
% Check if L1(1,1) and L1(1,2)
if L2(1,1) > 0 && L2(1,2) > 0
% If positive, break out of the loop
break;
end
end
%%
L1 = L1';
L2 = L2';
%%
XY1=[xCoordinates, yCoordinates, xCoordinates.^0];
G1=XY1*L1<0; %Group points into 3 regions.
G3=XY1*L2>0;
G2=~(G1|G3);
fn=@(G) mkpshape(xCoordinates, yCoordinates,G);
p(1)=fn(G1); %form polyshapes
p(2)=fn(G2);
p(3)=fn(G3);
%Visualize
h=plot(p); hold on
[h.FaceColor]=deal('g','b','m');
[h.EdgeColor]=deal('none');
plot(xCoordinates,yCoordinates,'.'); hold off
axis([1.263, 2.263 0.642, 1.442])
axis equal
view(90, -90); %rotation
Matt J
on 17 Jan 2024
You should check whether the points were properly separated by plotting
plot(xCoordinates(G1), yCoordinates(G1))
and similarly for G2 and G3.
Debora Baumann
on 17 Jan 2024
Matt J
on 17 Jan 2024
You're quite welcome, but if you consider your question resolved please Accept-click the answer.
Debora Baumann
on 18 Jan 2024
It's me again. It almost works now (even my imline), but i cannot find the reason why this is happening:
I plotted everything extra and the points seem to be sorted correctly. It's just these four points that mess everything up. Could you run the script below and see if it's happening for you too? I checked my line repeatedly, and it did not help :( do you find the issue?
clear all
%% Points outline
% Assuming you have a CSV file with coordinate points, load the data
coordinateData = readtable('outline_points_straight.csv'); % Replace with your actual file name
%ensure that name of stored G1, G2, and G3 is correct!!
xCoordinates = coordinateData.x;
coordinateData.x = xCoordinates + 0.363;
yCoordinates = coordinateData.y;
coordinateData.y = yCoordinates + 4.342;
% Define conditions for xCoordinates and yCoordinates
xCondition = (coordinateData.x > 2.263) | (coordinateData.x < 1.263);
yCondition = (coordinateData.y > 1.442) | (coordinateData.y < 0.642);
% Delete rows that satisfy the conditions
coordinateData(xCondition | yCondition, :) = [];
xCoordinates = coordinateData.x;
yCoordinates = coordinateData.y;
%%
% Plot the coordinate points as a scatter plot with equal axes
figure;
scatter(xCoordinates, yCoordinates, 0.2, 'filled');
%
% % Add labels and title
xlabel('X');
ylabel('Y');
title('Shifted Coordinate Points');
axis equal; % Ensure equal axis scaling
view(90, -90); % Rotation
xlim([1.263, 2.263]);
ylim([0.642, 1.442]);
% Use imline to draw the first line and store its position
%double click on the line after final position and then copy points or draw second line and then enter
% Draw the first line
h1 = imline;
set(findobj(h1, 'Type', 'Line'), 'LineWidth', 0.2);
wait(h1);
endpoints1 = getPosition(h1);
% Draw the second line
h2 = imline;
set(findobj(h2, 'Type', 'Line'), 'LineWidth', 0.2);
wait(h2);
endpoints2 = getPosition(h2);
disp(endpoints1);
disp(endpoints2);
P1l = [endpoints1(1, 1), endpoints1(1, 2),1];
x_left = [endpoints1(1, 1), endpoints1(2, 1)];
y_left = [endpoints1(1, 2), endpoints1(2, 2)];
coefficients_left = polyfit(x_left, y_left, 1);
slope_left = coefficients_left(1);
intercept_left = coefficients_left(2);
% Display the equation of the line
disp(['Equation of the left line: y = ', num2str(slope_left), ' * x + ', num2str(intercept_left)]);
%%
% Create a vector with values between 0 and 5 with steps of 0.1
y_values = 0.642:0.01:1.442;
% Iterate over each y value
for y_left = y_values
% Solve for x
x_left = (y_left - intercept_left) / (slope_left);
% Calculate Cx and Cy
P2l=[x_left,y_left,1];
L1= cross(P1l,P2l);
% Check if L1(1,1) and L1(1,2)
if L1(1,1) > 0 && L1(1,2) > 0
% If positive, break out of the loop
break;
end
end
%% right line
P1r = [endpoints2(1, 1), endpoints2(1, 2),1];
%% % Starting and endpoint of line
x_right = [endpoints2(1, 1), endpoints2(2, 1)];
y_right = [endpoints2(1, 2), endpoints2(2, 2)];
coefficients_right = polyfit(x_right, y_right, 1);
slope_right = coefficients_right(1);
intercept_right = coefficients_right(2);
disp(['Equation of the right line: y = ', num2str(slope_right), ' * x + ', num2str(intercept_right)]);
%%
% Iterate over each y value
for y_right = y_values
% Solve for x
x_right = (y_right - intercept_right) / (slope_right);
% Calculate Cx and Cy
P2r=[x_right,y_right,1];
L2= cross(P1r,P2r);
% Check if L1(1,1) and L1(1,2)
if L2(1,1) > 0 && L2(1,2) > 0
% If positive, break out of the loop
break;
end
end
%%
L1 = L1';
L2 = L2';
%%
XY1 = [xCoordinates, yCoordinates, xCoordinates.^0];
G1 = XY1 * L1 < 0;
hold on;
plot(xCoordinates(G1), yCoordinates(G1), 'o', 'MarkerSize', 1, 'MarkerFaceColor', 'b');
axis([1.263, 2.263 0.642, 1.442]); %xmin, xmax, ymin, and ymax
axis equal
view(90, -90); %rotation
%%
G3=XY1*L2<0;
hold on;
plot(xCoordinates(G3), yCoordinates(G3), 'o', 'MarkerSize', 1, 'MarkerFaceColor', 'b');
axis([1.263, 2.263 0.642, 1.442]);
axis equal
view(90, -90); %rotation
%%
G2 = ~G1 & ~G3;
%G2=~(G1|G3);
hold on;
plot(xCoordinates(G2), yCoordinates(G2), 'o', 'MarkerSize', 1, 'MarkerFaceColor', 'b');
axis([1.263, 2.263 0.642, 1.442]);
axis equal
view(90, -90); %rotation
%%
% Save G1, G2, G3 to CSV files
save('G1_straightoZt.mat', 'G1');
save('G2_straightoZt.mat', 'G2');
save('G3_straightoZt.mat', 'G3');
%%
p3 = mkpshape(xCoordinates, yCoordinates, G3);
%%
fn=@(G) mkpshape(xCoordinates, yCoordinates,G);
p(1)=fn(G1); %form polyshapes
p(2)=fn(G2);
p(3)=fn(G3);
%Visualize
h=plot(p); hold on
%[h.FaceAlpha]=deal(1);
[h.FaceColor]=deal('k','b','r');
[h.EdgeColor]=deal('none');
plot(xCoordinates,yCoordinates,'.'); hold off
axis([1.263, 2.263 0.642, 1.442])
axis equal
view(90, -90); %rotation
Matt J
on 18 Jan 2024
I won't be able to reproduce what you're seeing because the manual step of using imline will be different for me than for you. I recommend you just attach your L1 and L2 in a .mat file.
Debora Baumann
on 18 Jan 2024
Debora Baumann
on 18 Jan 2024
Here's an altenative implementation of mkpshape, which uses this FEX download,
It seems to work quite well.
function p=mkpshape(x,y,G)
%Trace the boundary of the given points and convert to polyshapes
x=x(G); y=y(G);
pth = tspsearch([x,y],5);
p=polyshape(x(pth),y(pth),'Simplify',0);
end
Also, you need to fix the computation of G3,
G3=XY1*L2 > 0;
Debora Baumann
on 20 Jan 2024
More Answers (0)
Categories
Find more on Graphics Object Properties in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
