Пересечение эллипса и прямой

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

Скачать исходник + исполняемый файл

Эта задача возникает при разнообразных квази-3D, изометрических отрисовках, построении графиков, диаграмм. Ниже примеры задач, где мне это понадобилось.

Конструктор муссовых десертов. GDI+. Изометрия.
Диаграмма направленности антенны. GDI+. Изометрия

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

Теория

Рис.1. Постановка задачи

Формулировка

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

Rx и Ry – полуоси эллипса, радиусы по X и Y. Аналог a и b в теории.

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

(1) Уравнение эллипса

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

(2) Уравнение прямой

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

Выразим Y из формулы (2):

(3) Уравнение для Y

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

(4)

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

(5)

Подставим вместо b выражение (Y0k x0) , и получим в верхней части второго слагаемого такую запись:

kx – kx0 = k(x –x0).

Таким образом, формула (5) вырождается в:

(6)

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

(7) Константа

Преобразуем формулу (6), учитывая (7), к следующему виду:

x2 – 2 x0 x + (x02 – ω) = 0

Решаем квадратное уравнение и получаем:

Помним, что:

  • Если D < 0, корней нет;
  • Если D = 0, один корень;
  • Если D > 0, два корня.

Находим x и x1 по формулам :

(8) Итоговые уравнения для x и x1

Практика

Напишем функцию, реализующую представленный выше расчет.

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 по используемым типам:

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;

Пример использования

Заранее прошу прощения, что ради демонстрации всего лишь того, что функция рабочая, понаписал дополнительно еще кучу вещей. Но мне показалось полезным и нужным реализовать следующие вещи:

  1. Приятный и функциональный интерфейс, показывающий работу алгоритма в максимально большем количестве аспектов.
  2. Показать рабочий пример работы с GDI+, т.к. сайт в значительной степени ориентирован на него.
  3. Показать рабочий пример представленного здесь GDI+ Canvas. Продемонстрировать одну из возможностей, а именно – используя стандартный GDI и TCanvas, рисовать между тем с возможностями GDI+.

Клик мышки по полю со схемой блокирует/разблокирует реакцию на курсор. Т.е. линия, пересекающая эллипс, «застывает». Слева три режима отрисовки. Галка AntiAlias снимает/устанавливает режим сглаживания линий.

Пример использования функции

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

Опции проекта

P.S.

Меня опять терзают смутные сомнения… У Шпака — магнитофон, у посла — медальон…

Бунша И.В.

Функцию можно использовать, скажем, для нахождения координат точки на линии эллипса, отстоящей на angle градусов. Пишем что-то типа следующего:

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 – результат, точка на эллипсе.

Неужели нет другого способа? Это ж могучая тригонометрия, должен же быть способ сразу посчитать координаты. Смущает искусственность способа – по углу найти точку линии, а потом найти ее пересечение.

Конечно же такой способ есть. Это ж могучая тригонометрия.

Оставить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *