Intersection of areas via polyshape

Hi everyone,

I need a help with the calculation of the intersection area via polyshape. I wrote a code that is designed to calculate gamma stability area. So, when I run the code and then there exists a point in the area that I cannot delete. Could you help to delete that? I wrote the code as below and I attach the output figure.

Thanks,

clear
clc
close all

wnArr = linspace(125,250,400);
zetaArr = linspace(0.1,0.9,400);

complexRoot1 = zeros(1,length(wnArr));
complexRoot2 = zeros(1,length(wnArr));
complexRoot3 = zeros(1,length(wnArr));
complexRoot4 = zeros(1,length(wnArr));


wnLow = 125; wnUp = 250;
zetaLow = 0.1; zetaUp = 0.9;

for k = 1:length(wnArr)
    complexRoot1(k) = -zetaArr(k)*wnLow + 1i*wnLow*sqrt(1-zetaArr(k)^2);
    complexRoot2(k) = -zetaArr(k)*wnUp + 1i*wnUp*sqrt(1-zetaArr(k)^2);
    complexRoot3(k) = -zetaLow*wnArr(k) + 1i*wnArr(k)*sqrt(1-zetaLow^2);
    complexRoot4(k) = -zetaUp*wnArr(k) + 1i*wnArr(k)*sqrt(1-zetaUp^2);
end

real1 = real(complexRoot1); imag1 = imag(complexRoot1);
real2 = real(complexRoot2); imag2 = imag(complexRoot2);
real3 = real(complexRoot3); imag3 = imag(complexRoot3);
real4 = real(complexRoot4); imag4 = imag(complexRoot4);

vertice1 = [transpose([real1 real2 real3 real4]) transpose([imag1 imag2 imag3 imag4])];
vertice2 = [transpose([real1 real2 real3 real4]) transpose([-imag1 -imag2 -imag3 -imag4])];

pgon1 = polyshape([0 0 1 1],[1 0 0 1], 'Simplify',false);
pgon1.Vertices = vertice1;

pgon2 = polyshape([0 0 1 1],[1 0 0 1], 'Simplify',false);
pgon2.Vertices = vertice2;

plot(simplify(simplify(pgon1)),'FaceColor', 'c')
hold on
plot(simplify(simplify(pgon2)),'FaceColor', 'c')
hold off

Answer 1

[Nipun]

Hi Sena,

I understand that you require help with the calculation of the intersection area via polyshape while deleting a specific set of points from the shape.

I recommend using the "subtract" function in MATLAB to delete a particular polyshape from a given polyshape.

Here, I am attaching the mended code which should solve your query:

clear
clc
close all

wnArr = linspace(125, 250, 400);
zetaArr = linspace(0.1, 0.9, 400);

complexRoot1 = zeros(1, length(wnArr));
complexRoot2 = zeros(1, length(wnArr));
complexRoot3 = zeros(1, length(wnArr));
complexRoot4 = zeros(1, length(wnArr));

wnLow = 125; wnUp = 250;
zetaLow = 0.1; zetaUp = 0.9;

for k = 1:length(wnArr)
    complexRoot1(k) = -zetaArr(k) * wnLow + 1i * wnLow * sqrt(1 - zetaArr(k)^2);
    complexRoot2(k) = -zetaArr(k) * wnUp + 1i * wnUp * sqrt(1 - zetaArr(k)^2);
    complexRoot3(k) = -zetaLow * wnArr(k) + 1i * wnArr(k) * sqrt(1 - zetaLow^2);
    complexRoot4(k) = -zetaUp * wnArr(k) + 1i * wnArr(k) * sqrt(1 - zetaUp^2);
end

real1 = real(complexRoot1); imag1 = imag(complexRoot1);
real2 = real(complexRoot2); imag2 = imag(complexRoot2);
real3 = real(complexRoot3); imag3 = imag(complexRoot3);
real4 = real(complexRoot4); imag4 = imag(complexRoot4);

vertice1 = [transpose([real1 real2 real3 real4]) transpose([imag1 imag2 imag3 imag4])];
vertice2 = [transpose([real1 real2 real3 real4]) transpose([-imag1 -imag2 -imag3 -imag4])];

pgon1 = polyshape([0 0 1 1], [1 0 0 1], 'Simplify', false);
pgon1.Vertices = vertice1;

pgon2 = polyshape([0 0 1 1], [1 0 0 1], 'Simplify', false);
pgon2.Vertices = vertice2;

%% Create a circular exclusion region
exclusionRadius = 100;
exclusionCenter = [0, 100];
theta = linspace(0, 2 * pi, 100);
exclusionX = exclusionCenter(1) + exclusionRadius * cos(theta);
exclusionY = exclusionCenter(2) + exclusionRadius * sin(theta);

exclusionRegion = polyshape([0 0 1 1], [1 0 0 1]);

exclusionRegion.Vertices = [transpose(exclusionX), transpose(exclusionY)];

%% Exclude points inside the circular region
pgon1 = subtract(pgon1, exclusionRegion);
pgon2 = subtract(pgon2, exclusionRegion);

%% Plot the results
plot(pgon1, 'FaceColor', 'c');
hold on
plot(pgon2, 'FaceColor', 'c');
plot(exclusionRegion, 'EdgeColor', 'b', 'FaceColor', 'none', 'LineStyle','--');
hold off

Here is the code output:

Link to "subtract" documentation: https://www.mathworks.com/help/matlab/ref/polyshape.subtract.html

Hope this helps.

Regards,

Nipun

Answer 2

[Matt J]

pgon1 = union(regions(polyshape(vertice1,'Simplify',1)));

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

pgon2 = union(regions(polyshape(vertice2,'Simplify',1)));

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

plot([pgon1,pgon2],'FaceColor', 'c')