Точнее будет «Пересечение эллипса и прямой, проходящей через его центр». Слишком длинное название не хотелось выводить в меню и оглавление. Статья уже немного устарела, более универсальный случай рассмотрен в статье Пересечение эллипса и произвольной прямой.
Скачать исходник + исполняемый файл
Эта задача возникает при разнообразных квази-3D, изометрических отрисовках, построении графиков, диаграмм. Ниже примеры задач, где мне это понадобилось.


Перейдем непосредственно к задаче. Т.к. сайт ориентирован на программирование, и в первую очередь на Delphi, поэтому исходим из того, что эллипс задается прямоугольником, в который он вписан.
Теория

Формулировка
Есть эллипс, вписанный в некий прямоугольник. Точка [X0;Y0] – координаты центра. Есть точка в пространстве [X2;Y2]. Необходимо найти координаты [X;Y] и [X1;Y1], которые являются точками пересечения эллипса и прямой, проходящей через точки [X0;Y0] и [X2;Y2].
Rx и Ry – полуоси эллипса, радиусы по X и Y. Аналог a и b в теории.
Формула эллипса для нашего случая:

Формула прямой для нашего случая:

Нам известны все координаты, кроме X, Y. Решаем систему уравнений (1),(2).
Выразим Y из формулы (2):

Преобразуем (3) к привычному виду линейной зависимости:

Подставим (4) в формулу (1):

Подставим вместо b выражение (Y0 – k x0) , и получим в верхней части второго слагаемого такую запись:
kx – kx0 = k(x –x0).
Таким образом, формула (5) вырождается в:

В формуле (6) заменим результат на константу:

Преобразуем формулу (6), учитывая (7), к следующему виду:
x2 – 2 x0 x + (x02 – ω) = 0
Решаем квадратное уравнение и получаем:
Помним, что:
- Если D < 0, корней нет;
- Если D = 0, один корень;
- Если D > 0, два корня.
Находим x и x1 по формулам :

Практика
Напишем функцию, реализующую представленный выше расчет.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
function CalcIntersecLineEllipse (ARect: TxRect; p2: TxPoint; APrec : Extended = 0.000001): TxRect; var rx, ry : double; // радиусы dx, dy : double; // разницы между заданными координатами a,b,c : double; // коэффициенты квадратного уравнения x, x1 : double; // X-координаты пересечения прямой y, y1 : double; // Y-координаты пересечения прямой p0 : TxPoint; // координата центра k : double; // см.описание - коэфф. k w : double; // см.описание - значение w d : double; // см.описание - дискриминант D begin //-- это означает, что решения нет ------------------- result.Left := -1000000; //-- радиусы по X и Y -------------------------------- rx := (ARect.Right - ARect.Left) / 2; ry := (ARect.Bottom - ARect.Top) / 2; //-- при нулевом радиусе решать что-либо бессмысленно if (rx = 0) or (ry = 0) then exit; //-- координаты центра ------------------------------- p0.X := ARect.Left + rx; p0.Y := ARect.Top + ry; //-- считаем разности координат ---------------------- dx := p2.x - p0.x; dy := p2.y - p0.y; //-- обрабатываем "особые" случаи -------------------- if abs(dx) < APrec then begin if dy < 0 then result := xRect (p0.x, p0.y-ry, p0.x, p0.y+ry) else result := xRect (p0.x, p0.y+ry, p0.x, p0.y-ry); exit; end; if abs(dy) < APrec then begin if dx < 0 then result := xRect (p0.x-rx, p0.y, p0.x+rx, p0.y) else result := xRect (p0.x+rx, p0.y, p0.x-rx, p0.y); exit; end; //-- считаем коэфф. K см.формулу (4) ----------------- k := dy/dx; //-- считаем W в формуле (6) ------------------------- w := (1/sqr(Rx) + sqr(k)/sqr(ry)); if w = 0 then exit; w := 1/ w; //-- коэффициенты квадратного уравнения -------------- a := 1; b := -2*p0.x; c := sqr(p0.x) - w; //-- ситаем дискриминант ----------------------------- d := sqr(b) - 4*a*c; if d = 0 then exit; // корней нет //-- считаем координаты X ---------------------------- d := sqrt(d); x := (-b + d) / (2*a); x1 := (-b - d) / (2*a); //-- подставляем X в формулу (4) считаем координаты Y y := (x - p0.x) * (dy/dx) + p0.y; y1 := (x1 - p0.x) * (dy/dx) + p0.y; //-- записываем результат ---------------------------- a := abs(p2.X - x); b := abs(p2.X - x1); if (a < b) then result := xRect(x,y,x1,y1) else result := xRect(x1,y1,x,y); end; |
В коде присутствуют какие-то непонятные типы. Чтобы функция скомпилировалась, надо в предложение Uses добавить xIPDraw. Либо добавить сверху следующий код.
Выдержка из xIPDraw по используемым типам:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
Type //-- вещественная точка ------------------------------ PxPoint = ^TxPoint; TxPoint = packed record X : single; Y : single; end; //-- вещественный прямоугольник ---------------------- PxRect = ^TxRect; TxRect = packed record case Integer of 0: (Left, Top, Right, Bottom: single); 1: (TopLeft, BottomRight: TxPoint); end; function xRect (ALeft, ATop, ARight, ABottom : double) : TxRect; begin result.Top := ATop; result.Left := ALeft; result.Right := ARight; result.Bottom := ABottom; end; function xRect (ARect : TRect) : TxRect; begin result := xRect (ARect.Left,ARect.Top,ARect.Right,ARect.Bottom); end; |
Пример использования
Заранее прошу прощения, что ради демонстрации всего лишь того, что функция рабочая, понаписал дополнительно еще кучу вещей. Но мне показалось полезным и нужным реализовать следующие вещи:
- Приятный и функциональный интерфейс, показывающий работу алгоритма в максимально большем количестве аспектов.
- Показать рабочий пример работы с GDI+, т.к. сайт в значительной степени ориентирован на него.
- Показать рабочий пример представленного здесь GDI+ Canvas. Продемонстрировать одну из возможностей, а именно – используя стандартный GDI и TCanvas, рисовать между тем с возможностями GDI+.
Клик мышки по полю со схемой блокирует/разблокирует реакцию на курсор. Т.е. линия, пересекающая эллипс, «застывает». Слева три режима отрисовки. Галка AntiAlias снимает/устанавливает режим сглаживания линий.

Если Вы скачали и настроили GDI+ Canvas, подкаталог GDIPCanvas в примере становится неактуальным. Имеет смысл удалить его и убрать в опциях проекта SearchPath запись «GDIPCanvas\»:

P.S.
Меня опять терзают смутные сомнения… У Шпака — магнитофон, у посла — медальон…
Бунша И.В.
Функцию можно использовать, скажем, для нахождения координат точки на линии эллипса, отстоящей на angle градусов. Пишем что-то типа следующего:
1 2 3 4 |
SinCos(angle * pi/180, sn,cs); p.X := cnt.X + round(w * cs); p.Y := cnt.Y + round(w * sn); v := CalcIntersecLineEllipse(rct,p); |
- sn, cs – синус и косинус угла angle (в градусах).
- Cnt – координаты центра
- W – ширина прямоугольника, описывающего эллипс
- P – точка линии, проходящей через центр.
- v – результат, точка на эллипсе.
Неужели нет другого способа? Это ж могучая тригонометрия, должен же быть способ сразу посчитать координаты. Смущает искусственность способа – по углу найти точку линии, а потом найти ее пересечение.
Конечно же такой способ есть. Это ж могучая тригонометрия.