Connecting points to polygon

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:

1. filling the space between the lines so the geometry is represented clearly.
2. 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!

Answer 1

[Matt J]

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.

[Debora Baumann]

I only showed a snippet of the code, this is the complete code:

clear all
clc
%storing files in ParaView
%file name: Slice...
%
%change following variables
slice_level = 0.7775;
simulation_nr = 4;
form = 'straight';
file_name = sprintf('slice%.4f_U_grid_V%d_%s.csv', slice_level, simulation_nr, form);
%read files and rename variable
folder_path = 'C:\Users\debor\Documents\Masterarbeit\Simulation\Versuch4_oZ\Slice0.7775';
files = dir(fullfile(folder_path,'*.csv'));

for i=1:length(files)
    file_path=fullfile(folder_path,files(i).name);
    variable_name = genvarname(['csv_' num2str(i)]);
    eval([variable_name '=readtable(file_path);']);

    prefix = num2str(i/10);
   eval([variable_name ' = renamevars(' variable_name ','...
        '{''Time'',''Points_0'',''Points_1'',''Points_2'',''U_0'',''U_1'',''U_2''},'...
        '{''' prefix ''',''x_' prefix ''',''z_' prefix ''',''y_' prefix ''',''Ux_' prefix ''',' ...
        '''Uz_' prefix ''',''Uy_' prefix '''});']);
end

%combine files into one file
allCSV = table();
for i=1:length(files) 
    variable_name = ['csv_' num2str(i)];
    allCSV = horzcat(allCSV, eval(variable_name));
end
%writetable(allCSV,'combi_allCSV.csv');
%% remove times und Points and shift left lower corner to coordinate system origin
%check if variables are accurate
U_table_shift = removevars(allCSV,["0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1", ...
    "1.1","1.2","x_0.2","y_0.2","z_0.2","x_0.3","y_0.3","z_0.3","x_0.4","y_0.4","z_0.4","x_0.5","y_0.5","z_0.5","x_0.6","y_0.6","z_0.6","x_0.7","y_0.7","z_0.7", ...
    "x_0.8","y_0.8","z_0.8","x_0.9","y_0.9","z_0.9","x_1","y_1","z_1","x_1.1","y_1.1","z_1.1","x_1.2","y_1.2","z_1.2"]);

Error using tabular/parenDelete
Unrecognized table variable name '0.1'.

Error in tabular/removevars (line 28)
a(:,vars) = [];

% U_table_shift = removevars(allCSV,["0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1", ...
%     "1.1","1.2","1.3","x_0.2","y_0.2","z_0.2","x_0.3","y_0.3","z_0.3","x_0.4","y_0.4","z_0.4","x_0.5","y_0.5","z_0.5","x_0.6","y_0.6","z_0.6","x_0.7","y_0.7","z_0.7", ...
%     "x_0.8","y_0.8","z_0.8","x_0.9","y_0.9","z_0.9","x_1","y_1","z_1","x_1.1","y_1.1","z_1.1","x_1.2","y_1.2","z_1.3","x_1.3","y_1.3"]);
%writetable (allCSVedited,'C:\Users\debor\Documents\Masterarbeit\Simulation\Versuch5_oZ\combiU_V5_slice.csv');

% Add 0.363 to each cell in the 'x_0.1' column
U_table_shift.("x_0.1") = U_table_shift.("x_0.1") + 0.363;

% Add 4.342 to each cell in the 'y_0.1' column
U_table_shift.("y_0.1") = U_table_shift.("y_0.1") + 4.342;
%% 

% Find columns starting with 'Ux', 'Uy', and 'Uz'
columns_start_Ux = startsWith(U_table_shift.Properties.VariableNames, 'Ux');
columns_Ux = find(columns_start_Ux);
columns_start_Uy = startsWith(U_table_shift.Properties.VariableNames, 'Uy');
columns_Uy = find(columns_start_Uy);
columns_start_Uz = startsWith(U_table_shift.Properties.VariableNames, 'Uz');
columns_Uz = find(columns_start_Uz);

%Find columns with point information
columns_start_x = startsWith(U_table_shift.Properties.VariableNames, 'x');
columns_points_x = find(columns_start_x);
columns_start_y = startsWith(U_table_shift.Properties.VariableNames, 'y');
columns_points_y = find(columns_start_y);
columns_start_z = startsWith(U_table_shift.Properties.VariableNames, 'z');
columns_points_z = find(columns_start_z);


%% 
% Plot 3D scatter plot
% figure;
% scatter(U_table_shift.("x_0.1"), U_table_shift.("y_0.1"), 'filled');
% title('Scatter Plot of Points in Slice');
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% grid on;
% axis equal;
% view (90,-90) %rotation
%% calculating mean 
%muss mit NMSE und Moving average bestimmt werden (ab wann NMSE klein und
%MA halbwegs konstant?)
prefix_condition = 0.6;
% Identify columns to keep based on the prefix condition
columns_to_keep = U_table_shift(:, str2double(extractAfter(U_table_shift.Properties.VariableNames, 3)) >= prefix_condition);
columns_to_keep = U_table_shift(:, str2double(extractAfter(U_table_shift.Properties.VariableNames, 3)) >= prefix_condition);
columns_to_keep = U_table_shift(:, str2double(extractAfter(U_table_shift.Properties.VariableNames, 3)) >= prefix_condition);
%% table with valid velocity values
columns_to_keep_x = U_table_shift(:, columns_points_x);
columns_to_keep_y = U_table_shift(:, columns_points_y);
columns_to_keep_z = U_table_shift(:, columns_points_z);

% Rename the columns to 'x', 'y', and 'z'
columns_to_keep_x.Properties.VariableNames(startsWith(columns_to_keep_x.Properties.VariableNames, 'x_0.1')) = {'x'};
columns_to_keep_y.Properties.VariableNames(startsWith(columns_to_keep_y.Properties.VariableNames, 'y_0.1')) = {'y'};
columns_to_keep_z.Properties.VariableNames(startsWith(columns_to_keep_z.Properties.VariableNames, 'z_0.1')) = {'z'};

U_table_valid = [columns_to_keep_x, columns_to_keep_y, columns_to_keep_z, columns_to_keep];
%% mean calculation
column_names1 = U_table_valid.Properties.VariableNames;

% Extract column indices from the table
columns_to_keep_Ux = find(startsWith(column_names1, 'Ux'));
columns_to_keep_Uy = find(startsWith(column_names1, 'Uy'));
columns_to_keep_Uz = find(startsWith(column_names1, 'Uz'));

% Calculate mean for each row based on selected columns, trans for
% transient
U_mean_trans_x = mean(U_table_valid{:, columns_to_keep_Ux}, 2);
U_mean_trans_y = mean(U_table_valid{:, columns_to_keep_Uy}, 2);
U_mean_trans_z = mean(U_table_valid{:, columns_to_keep_Uz}, 2);

% Add means to the table
U_table_valid.U_mean_trans_x = U_mean_trans_x;
U_table_valid.U_mean_trans_y = U_mean_trans_y;
U_table_valid.U_mean_trans_z = U_mean_trans_z;
%% 

%magnitude of velocity 
U_abs_length = sqrt(U_table_valid.U_mean_trans_x.^2 + ...
    U_table_valid.U_mean_trans_y.^2 + U_table_valid.U_mean_trans_z.^2);

%adding it to the table
U_table_valid.U_abs_length = U_abs_length;
%writetable(U_table_valid, file_name);
%% Visualization 
% Create a scatter plot of entire simulation space
figure;
scatter(U_table_valid.("x"), U_table_valid.("y"), 5, U_table_valid.U_abs_length, 'filled');

% Add labels and title
xlabel('X');
ylabel('Y');
%zlabel('Z');
title('3D Scatter Plot Colored by U\_abs\_length');
colorbar;  % Add colorbar to show U_abs_length values
view (90,-90) %rotation
axis equal;
hold on

%% Points outline
%% Points outline
% Assuming you have a CSV file with coordinate points, load the data
coordinateData = readtable('C:\Users\debor\Documents\Masterarbeit\Simulation\outline_points_straight.csv'); % Replace with your actual file name

% Extract x and y coordinates from the table
xCoordinates = coordinateData.x;
yCoordinates = coordinateData.y;

% Create a new table
geometry_points_outline = table(xCoordinates, yCoordinates);

% Shift the coordinates
geometry_points_outline.xCoordinates = geometry_points_outline.xCoordinates + 0.363;
geometry_points_outline.yCoordinates = geometry_points_outline.yCoordinates + 4.342;

% Delete rows based on x-Coordinates condition
xCondition = (geometry_points_outline.xCoordinates > 2.263) | (geometry_points_outline.xCoordinates < 1.263);
geometry_points_outline(xCondition, :) = [];

%Delete rows based on y-Coordinates condition
yCondition = (geometry_points_outline.yCoordinates > 1.442) | (geometry_points_outline.yCoordinates < 0.642);
geometry_points_outline(yCondition, :) = [];
writetable(geometry_points_outline, 'geometry_points_outline.csv');
% % Plot the coordinate points as a scatter plot with equal axes
% figure;
% scatter(geometry_points_outline.xCoordinates, geometry_points_outline.yCoordinates, 5, 'filled');
% 
% % Add labels and title
% xlabel('X');
% ylabel('Y');
% title('Shifted Coordinate Points');
% axis equal;  % Ensure equal axis scaling
% view (90,-90) %rotation
%% 

% Compute pairwise distances between points
distances = pdist2([geometry_points_outline.xCoordinates, geometry_points_outline.yCoordinates], [geometry_points_outline.xCoordinates, geometry_points_outline.yCoordinates]);

% Find the indices of the two nearest neighbors for each point
[~, nearestIndices] = mink(distances, 3, 2);  % Find the 3 nearest neighbors (including itself)

% Plot lines between the two closest points
figure;
hold on;
for i = 1:height(geometry_points_outline)
    nearestPoints = nearestIndices(i, :);
    
    % Plot lines between the two closest points in black
    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

% Add labels and title
xlabel('X');
ylabel('Y');
title('Lines Between Two Closest Points');
axis equal;
view (90,-90) %rotation
xlim([1.263, 2.263]);
ylim([0.642, 1.442]);



%%

% Define fine grid spacing
gridSpacing = 0.01;  % Specify the desired grid spacing

% Define the range of x and y based on your data
% minX = min(U_table_valid.("x"));
% maxX = max(U_table_valid.("x"));
% minY = min(U_table_valid.("y"));
% maxY = max(U_table_valid.("y"));

% Define the range of x and y based on wanted image section
minX =  0.9+0.363;
maxX =  1.9+0.363;
minY = -3.7+4.342;
maxY = -2.9+4.342;
% Create meshgrid with fixed grid spacing
[X, Y] = meshgrid(minX:gridSpacing:maxX, minY:gridSpacing:maxY);
% Store fine meshgrid coordinates in a table
fineGridPoints = table(X(:), Y(:), 'VariableNames', {'x_fine', 'y_fine'});
writetable(fineGridPoints, 'fineGridPoints.csv');
% Plot the grid points
%figure;
%scatter(X(:), Y(:), 5, 'k', 'filled');  % Scatter plot of grid points

% Add labels and title
% xlabel('X');
% ylabel('Y'); title('Grid Points');

% Define the coarser grid spacing for velocity vectors
coarseGridSpacing = 0.05;  % Specify the desired coarser grid spacing

% Create coarser meshgrid
[X_coarse, Y_coarse] = meshgrid(minX:coarseGridSpacing:maxX-coarseGridSpacing, minY+coarseGridSpacing:coarseGridSpacing:maxY-coarseGridSpacing);
% Store coarser meshgrid coordinates in a table
coarseGridPoints = table(X_coarse(:), Y_coarse(:), 'VariableNames', {'x_coarse', 'y_coarse'});
%% 
% Interpolate U_abs_length onto the fine grid
Ux_fine = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_mean_trans_x, X, Y);
Uy_fine = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_mean_trans_y, X, Y);
Uz_fine = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_mean_trans_z, X, Y);
Umag_fine = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_abs_length, X, Y);

% Interpolate U_abs_length, U_mean_trans_x, and U_mean_trans_y onto the coarser grid
Z_coarse = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_abs_length, X_coarse, Y_coarse);
Ux_coarse = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_mean_trans_x, X_coarse, Y_coarse);
Uy_coarse = griddata(U_table_valid.("x"), U_table_valid.("y"), U_table_valid.U_mean_trans_y, X_coarse, Y_coarse);

% Create a new table for fine grid interpolation
fineGridInterpolationTable = table(X(:), Y(:), Ux_fine(:), Uy_fine(:),Uz_fine(:), Umag_fine(:), 'VariableNames', {'x_fine', 'y_fine', 'Ux_fine', 'Uy_fine','Uz_fine', 'Umag_fine'});

% Construct the table name based on file_name
tableName = ['fineGrid_' file_name];

% Write the table to a CSV file
writetable(fineGridInterpolationTable, [tableName]);
%% 
figure;
hold on;

% Contour Plot of U
contourLevels = 0:0.05:max(U_table_valid.U_abs_length);

% Create contour plot of velocity magnitude
contourf(X, Y, Umag_fine, contourLevels);

% Plot lines between the two closest points
for i = 1:height(geometry_points_outline)
    nearestPoints = nearestIndices(i, :);

    % Plot lines between the two closest points in black
    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


% Add velocity vectors using quiver
quiver(X_coarse, Y_coarse, Ux_coarse, Uy_coarse, 'k', 'AutoScaleFactor', 2);  % Plot velocity vectors

% Add labels and title
xlabel('X');
ylabel('Y');
title('Combined Plot: Lines Between Two Closest Points and Contour Plot of U');
xlim([1.263, 2.263]);
ylim([0.642, 1.442]);

% Optionally, add colorbar for contour levels
colorbar('Ticks', contourLevels, 'TickLabels', num2str(contourLevels', '%.2f'));
view(90, -90);  %rotation
axis equal;
xticks(0:0.1:max(X(:)));
yticks(0:0.1:max(Y(:)));
colormap(jet);
hold off

%% 

I'm not sure how polyshape will work with my points, since they don't have he correct order.

[Matt J]

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]

Sorry, i hope it works now. And thanks for trying :)

plotgeometry()

Error using readtable
Unable to find or open '...outline_points_straight.csv'. Check the path and filename or file permissions.

Error in plotgeometry (line 3)
coordinateData = readtable('...outline_points_straight.csv'); % Replace with your actual file name

[Matt J]

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]

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]

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]

Thanks for your elaborate answer :) I'm still struggling (and feeling very incompetent) to get the different L1 and L2. I'm not getting the imline function to work properly for me. I can only draw one line and the endpoints are stored the second i finish my first attempt. Currently the code section looks like this (i stopped attempting to get more than one line after a while, therefore here is my basic approach with imline):

%% 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;

% Plot the coordinate points as a scatter plot with equal axes

figure;

scatter(geometry_points_outline.xCoordinates, geometry_points_outline.yCoordinates, 1, '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

h1 = imline;

endpoints1=position(h1);

%%

% Starting and endpoint from line

x = [endpoints1(1, 1), endpoints1(2, 1)];

y = [endpoints1(1, 2), endpoints1(2, 2)];

coefficients = polyfit(x, y, 1);

slope = coefficients(1);

intercept = coefficients(2);

% Display the equation of the line

disp(['Equation of the line: y = ', num2str(slope), ' * x + ', num2str(intercept)]);

I had to shift the x and y coordinates, because i have to fit them into another plot later. I don't really understand how you have 3 numbers in L1 and L2. Where does the third number come from, because from my understanding there is only the interception with the y-axis and the slope.

I hope your patience has not run out and you can help me to fully understand the code!

[Matt J]

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]

That makes a lot of sense! Now I'm receiving following error:

Index in position 1 is invalid. Array indices must be positive integers or logical values.

Error in mkpshape (line 5)
 V(end,:)=[];

Error in polygons_geometry (line 61)
fn=@(G) mkpshape(xCoordinates, yCoordinates,G);

Error in polygons_geometry (line 63)
p(2)=fn(G2);

My current code:

% 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.99159922178988 0.63927626459144;1.31455642023346 1.10814785992218];

endpoints_oZ_straight_right= [2.23533829024068 0.975363085155054;1.87795187943582 1.46243518189202];

%line follows equation A*x+B*y+C=0 --> can be represented as L*[x;y;1]=0

%L is the coeffcient vector L=[A,B,C]

% Starting and endpoint from line

% x = [endpoints1(1, 1), endpoints1(2, 1),1];

% y = [endpoints1(1, 2), endpoints1(2, 2),1];

x1 = [endpoints_oZ_straight_left(1, 1), endpoints_oZ_straight_left(2, 1),1];

y1 = [endpoints_oZ_straight_left(1, 2), endpoints_oZ_straight_left(2, 2),1];

L1 = cross(x1,y1);

L1 = L1';

x2 = [endpoints_oZ_straight_right(1, 1), endpoints_oZ_straight_right(2, 1),1];

y2 = [endpoints_oZ_straight_right(1, 2), endpoints_oZ_straight_right(2, 2),1];

L2 = cross (x2,y2);

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])

Is the issue that i have negative numbers in my L1 and L2?

[Matt J]

The only way to know exactly is to trap the error with dbstop and check what V contains.

However, I think you probably do want to scale L1 and L2 to have positive numbers. At least L1(1:2) and L2(1:2) should be positive.

[Debora Baumann]

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]

You should check whether the points were properly separated by plotting

plot(xCoordinates(G1), yCoordinates(G1))

and similarly for G2 and G3.

[Debora Baumann]

The issue was definitly my line, the imline is still not properly working for me, but i can copy the endpoints, so it's fine :) Thanks a lot for your help!

[Matt J]

You're quite welcome, but if you consider your question resolved please Accept-click the answer.

[Debora Baumann]

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]

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]

Right, here they are.

[Debora Baumann]

I think you have to add the following :) then they look normal

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;

[Matt J]

Here's an altenative implementation of mkpshape, which uses this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/71226-tspsearch

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]

Amazing, thanks for your effort!