Точки пересечения эллипса и прямой
Собственно цель в названии вопроса. У меня имеются координаты верхнего левого угла и нижнего правого угла ограничивающего прямоугольника для эллипса (другими словами эллипс вписан в этот прямоугольник). Также, имеются две точки прямой. Координатная система - левый верхний угол окна 0;0, y направлена вниз, x направлена вправо. Как пытался решить:
Вот уравнение эллипса с центром x0;y0 (мой случай):
Вот уравнение прямой:
Взял уравнение эллипса и прямой, упростил, выразил из них X. Через дискриминант получаю:
- если
D < 0- корней нет, прямая не пересекает эллипс; - если
D == 0- корень один, прямая является касательной к эллипсу (могу быть не прав); - если
D > 0- корней два, как раз точки пересечения прямой эллипса.
И еще один исключительный случай: если у прямой x2 - x1 очень близко к нулю, то формула меняется, и через дискриминант мы уже выражаем y, а x будет равен x1 || x2 (так как они совпадают, в моем случае координаты приходят целочисленные).
Я проверил этот способ, он отлично работает с эллипсом который является окружностью. Однако, как только у эллипса a != b - алгоритм ошибается, и я никак не могу найти проблему. Ниже прикладываю подробный код (специально расписал длиннее, чтобы было легче найти проблему), в котором заложен вышеописанный алгоритм + исключительной случай, когда x2 - x1 очень близко к нулю:
private void calculatePointsIntersectBlock(Rectangle rect, Point pt1, Point pt2){
Rectangle rectCoordinate = new Rectangle(rect.x, rect.y, rect.x+ rect.width, rect.y+ rect.height);// координаты эллипса
double _a = rect.width / 2.0; // полуось эллипса a
double _b = rect.height / 2.0; // полуось эллипса b
double x0 = (rectCoordinate.x + rectCoordinate.width) /2.0; // x центр эллипса
double y0 = (rectCoordinate.y + rectCoordinate.height) /2.0; // y центр эллипса
// уравнение прямой
double k = (pt2.y - pt1.y) / (double)(pt2.x - pt1.x);
double b = pt1.y - k * pt1.x;
// уравнение эллипса
double v = (_a*_a)*(_b*_b);
double s = b - y0;
// собранное квадратное уравнение
double A = -1;
double B = -1;
double C = -1;
if(pt2.x - pt1.x != 0) {
A = (_a * _a) + (_b * _b) * (k * k);
B = 2 * ((_a * _a) * k * s - (_b * _b) * x0);
C = (_a * _a) * (x0 * x0) + (_b * _b) * (s * s) - v;
}else{ // исключительная ситуация x2 - x1 близко к нулю
double w = (1 - (Math.pow(pt1.x - x0, 2) / (_a*_a))) * (_b*_b);
A = 1;
B = -2 * y0;
C = (y0 * y0) - w;
}
// решаем и выводим два корня
double D = (B * B) - 4*A*C;
// подставляем и находим x и y
if(pt2.x - pt1.x != 0) {
// одна точка пересечения
if (D < 0) {
return; // точек пересечения нет
} else if (D == 0) {
double x1 = (-B * 2) / (2 * A);
Point single = new Point((int) x1, (int) (k * x1 + b)); // единственная точка
} else {
// две точки
double x1 = (-B - Math.sqrt(D)) / (2 * A);
double x2 = (-B + Math.sqrt(D)) / (2 * A);
double y1 = k * x1 + b;
double y2 = k * x2 + b;
Point first = new Point((int) x1, (int) y1); // первая точка
Point second = new Point((int) x2, (int) y2); // вторая точка
}
}else{ // исключительная ситуация x2 - x1 близко к нулю
if (D < 0) {
return; // точек пересечения нет
} else if (D == 0) {
// одна точка
double y1 = (-B * 2) / (2 * A);
Point single = new Point((int) 0, (int) y1); // единственная точка
} else {
// две точки
double y1 = (-B - Math.sqrt(D)) / (2 * A);
double y2 = (-B + Math.sqrt(D)) / (2 * A);
double x1 = pt1.x;
double x2 = pt1.x;
Point first = new Point((int) x1, (int) y1); // первая точка
Point second = new Point((int) x2, (int) y2); // вторая точка
}
}
}
Помогите найти ошибку. Проверял на бумаге, но никак не могу найти, вроде бы всё сходится.
Ответы (1 шт):
Ошибки
Я буду использовать формат diff. Центр эллипса вычислен с ошибкой:
< double x0 = (rectCoordinate.x + rectCoordinate.width) /2.0; // x центр эллипса
< double y0 = (rectCoordinate.y + rectCoordinate.height) /2.0; // y центр эллипса
---
> double x0 = rectCoordinate.x + rectCoordinate.width / 2.0; // x центр эллипса
> double y0 = rectCoordinate.y + rectCoordinate.height / 2.0; // y центр эллипса
В коэффициентах A и C перепутаны полуоси эллипса:
< A = (_a * _a) + (_b * _b) * (k * k);
---
> A = (_b * _b) + (_a * _a) * (k * k);
< C = (_a * _a) * (x0 * x0) + (_b * _b) * (s * s) - v;
---
> C = (_b * _b) * (x0 * x0) + (_a * _a) * (s * s) - v;
После исправления ошибок программа вычисляет точки пересечения, но чтобы убедится в этом без использования отладчика, нужно устранить ...
Недостатки
Функция ничего не возвращает, проверить её можно только в отладчике. Это возможно, но долго и муторно. Исправим возвращаемое значение функции:
< private void calculatePointsIntersectBlock(Rectangle rect, Point pt1, Point pt2){
---
> public static DPoint[] calculatePointsIntersectBlock(Rectangle rect, Point pt1, Point pt2) {
DPoint - точка с координатами типа double. Типом Point я пользоваться не хочу, он целочисленный, а мне нужны неокруглённые значения для проверки.
Возвращаемые значения будут такими:
return new DPoint[] {}; // точек пересечения нет
...
return new DPoint[] {single}; // единственная точка
...
return new DPoint[] {first, second}; // пара точек
Теперь можно написать тесты (14 штук достаточно чтобы проверить все случаи) и убедится что исправленный код работает как требуется.
Комментарий
Ошибки из первого пункта возникли так как код сложный, а его обоснования перед глазами нет. В этом случае необходим длинный пошаговый комментарий который объясняет как уравнения относятся к параметрам метода, как уравнения разрешаются и как строится ответ.
Упрощение кода
Код сложный: перевод параметров функции в величины, относящиеся к задаче, построение квадратного уравнения, решение уравнения, восстановление точек пересечения. Всё нужно разбить на функции, которые можно протестировать отдельно.
Уравнение прямой
Уравнение прямой вида y = kx + b как правило бесполезно в программировании. Во-первых, оно не годится для вертикальных прямых. Во-вторых, коэффициент k вычисляется как отношение (y2 - y1) / (x2 - x1). Почему это не хорошо - ниже.
Введите параметр t:
x = x1 + t(x2 - x1) y = y1 + t(y2 - y1)
Можно проверить что система уравнений выше задаёт прямую проходящую через точки (x1, y1), (x2, y2). Деления нет, специальных случаев нет. Эту систему уравнений следует подставить в уравнение эллипса, решить относительно t, по его значениям восстановить точки на прямой.
Код
Интерфейс подогнан под теорию, решение квадратного уравнения вынесено в отдельный метод, в комментариях изложена теория:
public static DPoint[] intersectEllipseLine(
double x0, // центр эллипса (x)
double y0, // центр эллипса (y)
double a , // горизонтальная полуось эллипса
double b , // вертикальная полуось эллипса
double x1, // первая точка определяющая прямую (x)
double y1, // первая точка определяющая прямую (y)
double x2, // вторая точка определяющая прямую (x)
double y2 // вторая точка определяющая прямую (y)
) {
// параметрическое (от t) уравнение прямой
// x = (x2 - x1)t + x1
// y = (y2 - y1)t + y1
//
// уравнение эллипса
// (x - x0)^2 / a^2 + (y - y0)^2 / b^2 = 1
// <->
// b^2 (x - x0)^2 + a^2 (y - y0)^2 = a^2 b^2
//
// подставляем прямую в эллипс
// b^2 ((x2 - x1)t + x1 - x0)^2 + a^2 ((y2 - y1)t + y1 - y0)^2 = a^2 b^2
// <-> (x21 := x2 - x1, y21 := y2 - y1, x10 := x1 - x0, y10 := y1 - y0)
// b^2 (x21 t + x10)^2 + a^2 (y21 t + y10)^2 = a^2 b^2
// <->
// (b^2 x21^2 + a^2 y21^2) t^2 +
// 2(b^2 x21 x10 + aa^2 y21 y10) t +
// (b^2 x10^2 + a^2 y10^2 - a^2 b^2) = 0
final double x21 = x2 - x1;
final double y21 = y2 - y1;
final double x10 = x1 - x0;
final double y10 = y1 - y0;
final double a2 = a * a;
final double b2 = b * b;
final double t[] = solveQuadraticEquation(
b2 * x21 * x21 + a2 * y21 * y21,
2 * (b2 * x21 * x10 + a2 * y21 * y10),
b2 * x10 * x10 + a2 * y10 * y10 - a2 * b2
);
// восстанавливаем точки пересечения
final DPoint p[] = new DPoint[t.length];
for (int i = 0; i < t.length; ++i) {
p[i] = new DPoint(
t[i] * (x2 - x1) + x1,
t[i] * (y2 - y1) + y1
);
}
return p;
}