# Intersection condition between two ellipses

조회 수: 21 (최근 30일)
Sathiya S V 2024년 7월 4일
편집: Bruno Luong 2024년 7월 5일
I have two ellipse using the ellipse equation. I need to know whether the two ellipses are intersecting or not as a Boolean value. True if intersecting, and False if not intersecting. The usage of Solve function is not working in my matlab version. So please suggest me a numerical way to write this code. I have copied the code I am using for it.
ellipseTwo = '((x+0)^2)/2^2 + ((y-5)^2)/4^2 = 1'; ellipseOne = '((x+0)^2)/2^2 + (y-5)^2/4^2 = 1';
ezplot(ellipseOne, [-10, 10, -10, 10]);
hold on
ezplot(ellipseTwo, [-10, 10, -10, 10]);
##### 댓글 수: 6이전 댓글 4개 표시이전 댓글 4개 숨기기
Sathiya 2024년 7월 4일
@Bruno Luong Circle radius is not direction dependent (can be considered as a scalar line defining the radius), but in the case of ellipse, the radius(major axis or minor axis or inbetween) depends on the direction, so it is vector. Aproximating ellipse to a circle, wont be good idea. Correct me if I am wrong.
Bruno Luong 2024년 7월 4일
편집: Bruno Luong 2024년 7월 4일
I don't think you get my idea. There is no approximation if tha aspect ratio of two ellipse are identical. By changing parameter, it becomes circle equation with new coordinates. This is NOT an approximation, but exact transformation.
If the aspect ratio is not identical (or radom tilted ellipses the problem become minimizing uadratic function under disc contraints; and ca, be solved with traust region algorithm. I won't go there yet since you have not confimed the asumption of aspect ratio.

댓글을 달려면 로그인하십시오.

### 답변 (4개)

Torsten 2024년 7월 4일
편집: Torsten 2024년 7월 4일
I'm surprised that we get a division by zero by the below solution method if the aspect ratio of the two ellipses is equal.
But I think sigma18 and sigma19 have a common divisor. So the problem might vanish if you use
sol = simplify(solve(e1subs,y,'MaxDegree',4))
sol = solve(e1subs,y,'MaxDegree',4)
below.
If all solutions are complex, the two ellipses don't intersect. If not, you have intersections.
If you invest a little effort, use pencil and paper and deduce the coefficients of the polynomial e1subs of degree 4. Then you can use the numerical "roots" command to deduce the possible intersections (and thus the information whether the ellipses intersect or not).
syms x y
syms x1 x2 y1 y2 real
syms a1 b1 a2 b2 real positive
e1 = (x-x1)^2/a1^2 + (y-y1)^2/b1^2 == 1;
e2 = (x-x2)^2/a2^2 + (y-y2)^2/b2^2 == 1;
solx = solve(e1*a1^2-e2*a2^2,x)
solx =
e1subs = subs(e1,x,solx)
e1subs =
sol = solve(e1subs,y,'MaxDegree',4)
sol =
solsubs = subs(sol,[a1 a2 b1 b2 x1 x2 y1 y2],[4 4 2 2 0 -5 5 2])
Error using sym/subs (line 156)
Division by zero.
solsubs(1:4)
##### 댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Torsten 2024년 7월 4일
편집: Torsten 2024년 7월 4일
There are still special cases for which the below code doesn't work, e.g. if Y1 = Y2.
% E1: (x-X1)^2/A1^2 + (y-Y1)^2/B1^2 = 1
% E2: (x-X2)^2/A2^2 + (y-Y2)^2/B2^2 = 1
%
% Numerical values for the ellipse parameters
A1 = 4;
A2 = 4;
B1 = 2;
B2 = 2;
X1 = 0;
X2 = -5;
Y1 = 5;
Y2 = 2;
%
syms x y
syms x1 x2 y1 y2 real
syms a1 b1 a2 b2 real positive
e1 = (x-x1)^2/a1^2 + (y-y1)^2/b1^2 == 1;
e2 = (x-x2)^2/a2^2 + (y-y2)^2/b2^2 == 1;
soly = solve(e1*b1^2-e2*b2^2,y);
e1subs = subs(e1,y,soly);
poly_x = fliplr(coeffs(lhs(e1subs)-rhs(e1subs),x));
sol_x = roots(subs(poly_x,[a1 a2 b1 b2 x1 x2 y1 y2],[A1 A2 B1 B2 X1 X2 Y1 Y2]));
poly_y = fliplr(coeffs(lhs(e1)-rhs(e1),y));
count = 0;
for i = 1:numel(sol_x)
if abs(imag(sol_x(i))) < 1e-8
sol_y = roots(subs(poly_y,[a1 a2 b1 b2 x1 x2 y1 y2 x],[A1 A2 B1 B2 X1 X2 Y1 Y2 sol_x(i)]));
for j = 1:2
if abs(imag(sol_y(j))) < 1e-8 & abs(subs(lhs(e2),[a1 a2 b1 b2 x1 x2 y1 y2 x y],[A1 A2 B1 B2 X1 X2 Y1 Y2 sol_x(i) sol_y(j)])-1) < 1e-8
count = count + 1;
SOL_X(count) = sol_x(i);
SOL_Y(count) = sol_y(j);
end
end
end
end
SOL_X
SOL_X =
SOL_Y
SOL_Y =
if isempty(SOL_X)
disp('Ellipses don''t intersect')
else
disp('Ellipses intersect')
end
Ellipses intersect

댓글을 달려면 로그인하십시오.

Aquatris 2024년 7월 4일
편집: Aquatris 2024년 7월 4일
Here is one way where accuracy will depend on selecting the threshold correctly. You simply solve the equation for y given an x vector, and then find the distance between the elips 1 points and elips 2 points. Then find if any distance is less than a threshold
x = -10:0.01:10;
% solve the equation for y and get index of imaginary solutions
y1 = sqrt(-((x+5).^2/4^2-1)*2^2);idx1 = imag(y1) ==0;
y2 = sqrt(-((x+0).^2/4^2-1)*2^2);idx2 = imag(y2) ==0;
% form x-y points of the elipses
x1 = x(idx1) ;x1 = [x1 flip(x1)];
y1 = y1(idx1);y1 = [y1 -flip(y1)]+2;
x2 = x(idx2) ;x2 = [x2 flip(x2)];
y2 = y2(idx2);y2 = [y2 -flip(y2)]+5;
% loopt over points and calculate their distance, this might be vectorized
for i = 1:length(x1)
for j = 1:length(x2)
dist(i,j) = sqrt((x1(i)-x2(j))^2+(y1(i)-y2(j))^2);
end
end
% find indeces where distance is below a threshold
[a,b] = find(dist<2e-3);
% plot
figure(2)
plot(x1,y1,x2,y2,x1(a),y1(a),'mx',x2(b),y2(b),'ko')
axis([-10 10 -10 10])
##### 댓글 수: 0이전 댓글 -2개 표시이전 댓글 -2개 숨기기

댓글을 달려면 로그인하십시오.

Bruno Luong 2024년 7월 4일
편집: Bruno Luong 2024년 7월 5일
Here is a solution for generic ellipses using FEX file here
% Ti and ci define two ellipses E1 and E2 of the form
% (x,y) st norm(Ti*Xi) <= 1 where Xi := [x;y]-ci
c1 = [0; 5]; % circle center #1
c2 = [-5; 2]; % circle center #2
T1 = [1/4 0;
0 1/2]; % Change variable matrix
T2 = T1;
T = T2/T1;
A = 2*(T'*T);
dc = c1-c2;
dc2 = T2*dc;
b = 2*T*dc2;
c = dc2'*dc2;
% TrustRegion.m File from FEX
% https://www.mathworks.com/matlabcentral/fileexchange/27596-least-square-with-2-norm-constraint
Xx = TrustRegion(A, b, 1);
% X = T1\Xx + c1
f = 1/2*Xx'*A*Xx+Xx'*b+c;
isintersect = f <= 1
##### 댓글 수: 0이전 댓글 -2개 표시이전 댓글 -2개 숨기기

댓글을 달려면 로그인하십시오.

Bruno Luong 2024년 7월 5일
편집: Bruno Luong 2024년 7월 5일
Use FSOLVE to find numerical solution. Only xaimum one can be found
c1 = [0; 5]; % circle center #1
c2 = [-5; 2]; % circle center #2
T1 = [1/4 0;
0 1/2]; % Change variable matrix
T2 = T1;
[X, fval, exitflag] = fsolve(@(X) myfun(X,c1,T1,c2,T2), randn(2,1))
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
X = 2x1
-1.8347 3.2228
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
fval = 2x1
1.0e-07 * 0.2337 0.2337
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
exitflag = 1
isintersect = norm(fval, Inf) < 1e-6 && exitflag >= 0
isintersect = logical
1
figure
hold on
h1=fimplicit(@(x,y) sum((T1*([x(:),y(:)]'-c1)).^2,1)-1,[-10 10]);
h2=fimplicit(@(x,y) sum((T2*([x(:),y(:)]'-c2)).^2,1)-1,[-10 10]);
h1.Color = 'k';
h2.Color = 'k';
plot(X(1),X(2),'ro','MarkerSize',10)
function F = myfun(X,c1,T1,c2,T2)
F = [sum((T1*(X-c1)).^2,1)-1;
sum((T2*(X-c2)).^2,1)-1];
end

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Mathematics에 대해 자세히 알아보기

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by