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

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

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

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

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

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

Перейдем непосредственно к задаче. Т.к. сайт ориентирован на программирование, и в первую очередь на 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 выражение (Y₀ – k x₀) , и получим в верхней части второго слагаемого такую запись:

kx – kx₀= k(x –x₀).

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

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

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

x² – 2 x₀ x + (x₀²– ω) = 0

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

Помним, что:

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

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

Практика

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

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;

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;

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 градусов. Пишем что-то типа следующего:

SinCos(angle * pi/180, sn,cs); 
p.X := cnt.X + round(w * cs); 
p.Y := cnt.Y + round(w * sn); 
v := CalcIntersecLineEllipse(rct,p);

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

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

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

0 0 голоса

Рейтинг статьи

0 комментариев

Межтекстовые Отзывы

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