I am a 9th grade student. In my spare time I do programming and one of my projects is my game. Thinking about the next update, I thought about how to improve the hitbox collision calculation system - almost all objects in the game had an oval/round shape, and describing rectangles around them, there was a noticeable error (especially when hitbox corners collided). So I decided to implement an oval hitbox
The task was to implement two functions:
- collide_ellipse_rect (Checks if an ellipse collides with a rectangle)
- collide_ellipse_ellipse (Checks if an ellipse collides with an ellipse)
I dealt with the first easily by substituting the four coordinates of the rectangle into the ellipse equation and checking the range of the resulting coordinate.
But with the second function, everything was not so great. I was looking for a solution that works quickly, in O(1), while with minimal error (because if the error was large, it would be easier to leave a rectangular hitbox). At first I used a geometric approach which was as follows:
- Calculate the radii of ellipses at the current angle
- Compare the distance between the centers of the ellpses and the sum of their radii:
- if the sum is greater than the distance, they intersect
- otherwise they do not
Unfortunately, this method does not work for prolate ellipses, but only for those tending to a circle.
Even then I understood that this problem has an algebraic solution, and since the maximum intersection of two ellipses is four, an equation of the fourth degree should appear. But due to a lack of knowledge in algebra, I could not bring the canonical equation of two ellipses to a polynomial of the fourth degree
After a very long search on the Internet for information about this, I still found a needle in a haystack. The wonderful user Mishamp left for me just the solution to this problem back in 2011
The original answer
при решенье пересичение двух елипсов випиши уравненние и реши их
если ти правильно все зделал то у тебя должно получится полином 4й степени и ето уравненние легко решаетса формулой Ферарри
я уже демонстрировал етот пример
(r1x,r1y) - радиус первого еллипса
(cx,cy) - координати второго еллипса
(r2x,r2y) - радиус второго еллипса
// для упрощение подсчетов виделем в памяти ячейку для квадратов, чтоби потом процесор не пересчитовал каждий раз const double cx_2 = cx * cx; const double cy_2 = cy * cy; const double r1x_2 = r1x*r1x; const double r1y_2 = r1y*r1y; const double r2x_2 = r2x*r2x; const double r2y_2 = r2y*r2y; const double RyRx = r2y_2/r2x_2 - r1y_2/r1x_2 ; const double z4 = -RyRx*RyRx; const double z3 = 4.0*cx*r2y_2*RyRx/r2x_2; const double z2 = - 4.0 * ( cy_2 * r2x_2 * r2y_2 + cx_2 * r2y_2 * r2y_2 )/(r2x_2 * r2x_2) - 2.0 * RyRx * (( cx_2 / r2x_2 - 1.0) * r2y_2 - cy_2 + r1y_2 ); const double z1 = (4.0 * cx * r2y_2 * (cy_2 * r2x_2 + r1y_2 * r2x_2 + (cx_2 - r2x_2) * r2y_2))/(r2x_2*r2x_2); const double z0 = 4.0 * cy_2 * r2y_2 - (4.0 * cx_2 * cy_2 * r2y_2)/r2x_2 - SQR(r1y_2 - cy_2 + (cx_2/r2x_2 - 1.0) * r2y_2);// решить уравненние: z4 * x^4 + z3 * x^3 + z2 * x^2 + z1 * x + z0 == 0
получем координати X пересичения
но чтоби получить координати Y пересечения нужно подставить X в уравненние первого елипса получем две точки но не обезательно обе являютса пересечениями поетому нужно координати X Y подставить в уравненния второго елипса если не подошли тогда ето не пересичение
Внимание!!! обратите внимание что могут пересичения пропадать когда они должни бить, ето иизза того что в формуле процесор вичислял сначала очень большое число больше чем 1e+10 (поетому округлил все слеледуюсчи цифри в ноль, помня толька старшии цифри ) потом делил и получалось число довольно приближенное к настоясчиму числу но ето не точнее число когда происходит проверка, точки отбрасиваютса иза етои неточности... формула не совпадает
как с етим боротса я даже не думал пока что.... мне небило необходимости но думаю длиная математика и проблеми не будет, хотя уверен можно обойтись условием поиска одной из двух точек подходясчий
Thus, we can write the algorithm in verbal form:
- Get the coefficients of the equation of the 4th degree
- Solve the equation and get no more than 4 X coordinates (but if the intersection ellipses coincide, there will be an infinite number)
- Substitute all X coordinates into the equation of the current ellipse and get 2 (or 1) Y coordinates. Thus, we got no more than 8 intersection points
- Substitute all the points in the equation of another ellipse, and if it is true, then this point is the intersection
Do not forget the case if one ellipse lies inside another, therefore, before calculations, we must check if the center of one of the ellipses lies inside another ellipse, then they definitely collide
Now I could write code, but also, I need to write functions for solving equations of all degrees from 1st to 4th. Unfortunately, we have not yet gone through the solution of equation 3rd and, especially, the 4th degree for the general case, so after a little thought, I decided to find ready-made algorithms. After a short search, I decided to take the algorithm from the online equation calculator. It was written in JavaScript and I could easily rewrite it in Python.
3rd degree equation calculator algorithm
Site from which the algorithm was taken. JS source:
function cubicsolve(dataForm)
{
var a = parseFloat(dataForm.aIn.value);
var b = parseFloat(dataForm.bIn.value);
var c = parseFloat(dataForm.cIn.value);
var d = parseFloat(dataForm.dIn.value);
if (a == 0)
{
alert("??????????? ???? ????? 0.");
return;
} //End if a == 0
if (d == 0)
{
alert("???? ?????? ????? 0.");
return;
} //End if d == 0
b /= a;
c /= a;
d /= a;
var disc, q, r, dum1, s, t, term1, r13;
q = (3.0*c - (b*b))/9.0;
r = -(27.0*d) + b*(9.0*c - 2.0*(b*b));
r /= 54.0;
disc = q*q*q + r*r;
dataForm.x1Im.value = 0; //The first root is always real.
term1 = (b/3.0);
if (disc > 0) { // one root real, two are complex
s = r + Math.sqrt(disc);
s = ((s < 0) ? -Math.pow(-s, (1.0/3.0)) : Math.pow(s, (1.0/3.0)));
t = r - Math.sqrt(disc);
t = ((t < 0) ? -Math.pow(-t, (1.0/3.0)) : Math.pow(t, (1.0/3.0)));
dataForm.x1Re.value = -term1 + s + t;
term1 += (s + t)/2.0;
dataForm.x3Re.value = dataForm.x2Re.value = -term1;
term1 = Math.sqrt(3.0)*(-t + s)/2;
dataForm.x2Im.value = term1;
dataForm.x3Im.value = -term1;
return;
}
// End if (disc > 0)
// The remaining options are all real
dataForm.x3Im.value = dataForm.x2Im.value = 0;
if (disc == 0){ // All roots real, at least two are equal.
r13 = ((r < 0) ? -Math.pow(-r,(1.0/3.0)) : Math.pow(r,(1.0/3.0)));
dataForm.x1Re.value = -term1 + 2.0*r13;
dataForm.x3Re.value = dataForm.x2Re.value = -(r13 + term1);
return;
} // End if (disc == 0)
// Only option left is that all roots are real and unequal (to get here, q < 0)
q = -q;
dum1 = q*q*q;
dum1 = Math.acos(r/Math.sqrt(dum1));
r13 = 2.0*Math.sqrt(q);
dataForm.x1Re.value = -term1 + r13*Math.cos(dum1/3.0);
dataForm.x2Re.value = -term1 + r13*Math.cos((dum1 + 2.0*Math.PI)/3.0);
dataForm.x3Re.value = -term1 + r13*Math.cos((dum1 + 4.0*Math.PI)/3.0);
return;
} //End of cubicSolve4th degree equation calculator algorithm
Site from which the algorithm was taken. JS source:
function quad4solve(dataForm)
{
var a = parseFloat(dataForm.aIn.value);
var b = parseFloat(dataForm.bIn.value);
var c = parseFloat(dataForm.cIn.value);
var d = parseFloat(dataForm.dIn.value);
var e = parseFloat(dataForm.eIn.value);
if (a == 0)
{
alert("The coefficient of the power four of x is 0. Please use the utility for a third degree quadratic.");
return;
}
if (e == 0)
{
alert("One root is 0. Now divide through by x and use the utility for a third degree quadratic to solve the resulting equation for the other three roots.");
return;
}
if (a != 1)
{
b /= a;
c /= a;
d /= a;
e /= a;
}
// Coefficients for cubic solver
var cb, cc, cd;
var discrim, q, r, RRe, RIm, DRe, DIm, dum1, ERe, EIm, s, t, term1, r13, sqR, y1, z1Re, z1Im, z2Re;
cb = -c;
cc = -4.0*e + d*b;
cd = -(b*b*e + d*d) + 4.0*c*e;
if (cd == 0)
{
alert("cd = 0.");
}
q = (3.0*cc - (cb*cb))/9.0;
r = -(27.0*cd) + cb*(9.0*cc - 2.0*(cb*cb));
r /= 54.0;
discrim = q*q*q + r*r;
term1 = (cb/3.0);
if (discrim > 0)
{
// 1 root real, 2 are complex
s = r + Math.sqrt(discrim);
s = ((s < 0) ? -Math.pow(-s, (1.0/3.0)) : Math.pow(s, (1.0/3.0)));
t = r - Math.sqrt(discrim);
t = ((t < 0) ? -Math.pow(-t, (1.0/3.0)) : Math.pow(t, (1.0/3.0)));
y1 = -term1 + s + t;
}
else
{
if (discrim == 0)
{
r13 = ((r < 0) ? -Math.pow(-r,(1.0/3.0)) : Math.pow(r,(1.0/3.0)));
y1 = -term1 + 2.0*r13;
}
else
{
q = -q;
dum1 = q*q*q;
dum1 = Math.acos(r/Math.sqrt(dum1));
r13 = 2.0*Math.sqrt(q);
y1 = -term1 + r13*Math.cos(dum1/3.0);
}
}
// Determined y1, a real root of the resolvent cubic.
term1 = b/4.0;
sqR = -c + term1*b + y1;
RRe = RIm = DRe = DIm = ERe = EIm = z1Re = z1Im = z2Re = 0;
if (sqR >= 0)
{
if (sqR == 0)
{
dum1 = -(4.0*e) + y1*y1;
if (dum1 < 0) //D and E will be complex
z1Im = 2.0*Math.sqrt(-dum1);
else
{ //else (dum1 >= 0)
z1Re = 2.0*Math.sqrt(dum1);
z2Re = -z1Re;
}
}
else
{
RRe = Math.sqrt(sqR);
z1Re = -(8.0*d + b*b*b)/4.0 + b*c;
z1Re /= RRe;
z2Re = -z1Re;
}
}
else
{
RIm = Math.sqrt(-sqR);
z1Im = -(8.0*d + b*b*b)/4.0 + b*c;
z1Im /= RIm;
z1Im = -z1Im;
}
z1Re += -(2.0*c + sqR) + 3.0*b*term1;
z2Re += -(2.0*c + sqR) + 3.0*b*term1;
//At this point, z1 and z2 should be the terms under the square root for D and E
if (z1Im == 0)
{ // Both z1 and z2 real
if (z1Re >= 0)
{
DRe = Math.sqrt(z1Re);
}
else
{
DIm = Math.sqrt(-z1Re);
}
if (z2Re >= 0)
{
ERe = Math.sqrt(z2Re);
}
else
{
EIm = Math.sqrt(-z2Re);
}
}
else
{
r = Math.sqrt(z1Re*z1Re + z1Im*z1Im);
r = Math.sqrt(r);
dum1 = Math.atan2(z1Im, z1Re);
dum1 /= 2; //Divide this angle by 2
ERe = DRe = r*Math.cos(dum1);
DIm = r*Math.sin(dum1);
EIm = -DIm;
}
dataForm.x1Re.value = -term1 + (RRe + DRe)/2;
dataForm.x1Im.value = (RIm + DIm)/2;
dataForm.x2Re.value = -(term1 + DRe/2) + RRe/2;
dataForm.x2Im.value = (-DIm + RIm)/2;
dataForm.x3Re.value = -(term1 + RRe/2) + ERe/2;
dataForm.x3Im.value = (-RIm + EIm)/2;
dataForm.x4Re.value = -(term1 + (RRe + ERe)/2);
dataForm.x4Im.value = -(RIm + EIm)/2;
return;
} As a result, I got the following code:
from math import sqrt, acos, cos, sin, atan2, pi
# Ellipse function with X argument
def fy(ellipse, x):
cx, cy, a, b = ellipse.__data__
d = 1 - (x - cx) ** 2 / a ** 2
if d >= 0:
d = sqrt(d)
return (cy - b * d, cy + b * d)
return ()
# Cubic root funciton
def cbrt(x):
return pow(x, 1 / 3) if x > 0 else -pow(-x, 1 / 3)
""" Start of equation solution functions """
def solve_p1(a, b):
if a == 0:
return set([0])
return set([-b / a])
def solve_p2(a, b, c):
# Linear equation
if a == 0:
return solve_p1(b, c)
D = b ** 2 - 4 * a * c
if D >= 0:
d = sqrt(D)
return set([
(-b - d) / 2 / a,
(-b + d) / 2 / a
])
return set()
def solve_p3(a, b, c, d):
# Quadratic equation
if a == 0:
return solve_p2(b, c, d)
if d == 0:
return set([
*solve_p2(a, b, c),
0
])
b /= a
c /= a
d /= a
q = (3 * c - b * b) / 9
r = (b * (9 * c - 2 * b * b) - 27 * d) / 54
D = q * q * q + r * r
term = b / 3
# One root is real (two are complex)
if D > 0:
s = cbrt(r + sqrt(D))
t = cbrt(r - sqrt(D))
return set([
s + t - term
])
# All roots are real (at least two are equal)
if D == 0:
r13 = cbrt(r)
return set([
2 * r13 - term,
-(r13 + term)
])
# All roots are real (all roots are unequal)
q = -q
dum1 = acos(r / sqrt(q * q * q))
r13 = 2 * sqrt(q)
return set([
r13 * cos(dum1 / 3) - term,
r13 * cos((dum1 + 2 * pi) / 3) - term,
r13 * cos((dum1 + 4 * pi) / 3) - term
])
def solve_p4(a, b, c, d, e):
# Cubic equation
if a == 0:
return solve_p3(b, c, d, e)
if e == 0:
return set([
*solve_p3(a, b, c, d),
0
])
b /= a
c /= a
d /= a
e /= a
cb = -c
cc = -4 * e + d * b
cd = -(b * b * e + d * d) + 4 * c * e
q = (3 * cc - cb * cb) / 9
r = (cb * (9 * cc - 2 * cb * cb) - 27 * cd) / 54
D = q * q * q + r * r
term = cb / 3
if D > 0:
s = cbrt(r + sqrt(D))
t = cbrt(r - sqrt(D))
y1 = s + t - term
else:
if D == 0:
r13 = cbrt(r)
y1 = 2 * r13 - term
else:
q = -q
dum1 = acos(r / sqrt(q * q * q))
r13 = 2 * sqrt(q)
y1 = r13 * cos(dum1 / 3) - term
term = b / 4
sqR = term * b + y1 - c
RRe = RIm = DRe = DIm = ERe = EIm = z1Re = z1Im = z2Re = 0
if sqR >= 0:
if sqR == 0:
dum1 = y1 * y1 - 4 * e
if dum1 < 0:
z1Im = 2 * sqrt(-dum1)
else:
z1Re = 2 * sqrt(dum1)
z2Re = - z1Re
else:
RRe = sqrt(sqR)
z1Re = -(8 * d + b * b * b) / 4 + b * c
z1Re /= RRe
z2Re = -z1Re
else:
RIm = sqrt(-sqR)
z1Im = -(8 * d + b * b * b) / 4 + b * c
z1Im /= RIm
z1Im = -z1Im
z1Re += 3 * b * term - 2 * c - sqR
z2Re += 3 * b * term - 2 * c - sqR
if z1Im == 0:
if z1Re >= 0:
DRe = sqrt(z1Re)
else:
DIm = sqrt(-z1Re)
if z2Re >= 0:
ERe = sqrt(z2Re)
else:
EIm = sqrt(-z2Re)
else:
r = sqrt(sqrt(z1Re * z1Re + z1Im * z1Im))
dum1 = atan2(z1Im, z1Re) / 2
ERe = DRe = r * cos(dum1)
DIm = r * sin(dum1)
EIm = -DIm
roots = (
((RRe + DRe) / 2 - term, (RIm + DIm) / 2),
(RRe / 2 - term - DRe / 2, (RIm - DIm) / 2),
(ERe / 2 - term - RRe / 2, (EIm - RIm) / 2),
(-term - (RRe + ERe) / 2, -(RIm + EIm) / 2)
)
# Filter real roots
return set([
i[0] for i in roots if abs(i[1]) < 1e-8
])
""" End of equation solution functions """
class Ellipse:
def __init__(self, x, y, a, b):
self.center = (x, y)
self.__data__ = [x, y, a, b]
def collide_point(self, x, y):
cx, cy, a, b = self.__data__
if (x - cx) ** 2 / a ** 2 + (y - cy) ** 2 / b ** 2 - 1e-6 <= 1:
return True
return False
def collide_ellipse(self, ellipse):
x1, y1, a1, b1 = self.__data__
x2, y2, a2, b2 = ellipse.__data__
if self.collide_point(*ellipse.center) or ellipse.collide_point(*self.center):
return True
# Transform the coordinates so that one of the ellipses is at the beginning
x2 -= x1
y2 -= y1
# Calculate the coefficients of the equation
x2_2 = x2 * x2
y2_2 = y2 * y2
a1_2 = a1 * a1
b1_2 = b1 * b1
a2_2 = a2 * a2
b2_2 = b2 * b2
RyRx = b2_2 / a2_2 - b1_2 / a1_2
coefs = (
-RyRx * RyRx,
4.0 * x2 * b2_2 * RyRx / a2_2,
-4.0 * (y2_2 * a2_2 * b2_2 + x2_2 * b2_2 * b2_2) / (a2_2 * a2_2) - 2.0 * RyRx * ((x2_2 / a2_2 - 1.0) * b2_2 - y2_2 + b1_2),
(4.0 * x2 * b2_2 * (y2_2 * a2_2 + b1_2 * a2_2 + (x2_2 - a2_2) * b2_2)) / (a2_2 * a2_2),
4.0 * y2_2 * b2_2 - (4.0 * x2_2 * y2_2 * b2_2) / a2_2 - (b1_2 - y2_2 + (x2_2 / a2_2 - 1.0) * b2_2)**2
)
# Solution of the equation and checking intersection points
for x in solve_p4(*coefs):
x += x1
for y in fy(self, x):
if ellipse.collide_point(x, y):
return True
return False
def test():
assert Ellipse(79, 73, 62, 45).collide_ellipse(Ellipse(87.5, 75, 45.5, 65))
assert not Ellipse(70, 45, 19, 35).collide_ellipse(Ellipse(145, 84.5, 60, 19.5))
assert not Ellipse(164, 291, 60, 90).collide_ellipse(Ellipse(92, 195.5, 42, 19.5))
assert Ellipse(75.5, 136, 55.5, 49).collide_ellipse(Ellipse(104, 122.5, 18, 22.5))
assert Ellipse(2.5, 44, 42.5, 12).collide_ellipse(Ellipse(59, 44, 15, 27))
assert Ellipse(157, 157, 99, 83).collide_ellipse(Ellipse(103, 116, 27, 25))
assert Ellipse(-42.5, -58, 57.5, 42).collide_ellipse(Ellipse(-42.5, -58, 57.5, 42))
assert Ellipse(39, 130, 5, 50).collide_ellipse(Ellipse(84, 87, 50, 5))
assert not Ellipse(45, 242, 40, 21).collide_ellipse(Ellipse(208, 242, 33, 27))
assert Ellipse(91, 152, 38, 60).collide_ellipse(Ellipse(146, 152, 17, 60))
# Test Ellipse
test()
# x, y - ellipse centers
# a, b - semiaxes of ellipses
x1, y1, a1, b1 = map(float, input().split())
x2, y2, a2, b2 = map(float, input().split())
e1 = Ellipse(x1, y1, a1, b1)
e2 = Ellipse(x2, y2, a2, b2)
print(e1.collide_ellipse(e2))The error is essentially large, but acceptable for my task, because one coordinate corresponds to one pixel. This solution works quickly, in O(1), although it is a big unit (however, this did not affect the current (60) fps in any way)
As it turned out, checking the collision of two ellipses is not such a trivial task, but I was able to find a solution to this problem and I am sharing it with you. I believe that this article can be useful not only in game development, but also in other areas. I hope it will help you to reduce the time spent on solving this problem.
YariKartoshe4ka, With love 08.02.2022