Точки пересечения кривых Безье

Точки пересечения кривых Безье

-
-
В статье описывается способ вычисления точек пересечения кривых Безье. Дело в том, что, например, для получения точек пересечения кривой Безье с прямой достаточно решить не очень сложное уравнение третьей степени. В случае двух кривых Безье получаем систему, которую, на первый взгляд, решить невозможно. Единственное, что удалось найти по этому вопросу в сети - исходник на JavaScript[1] с готовым решением (автор Kevin Lindsey). Мне же хотелось разобраться самому как получить решение. Хочу поделиться тем, что у меня получилось. Возможно, кому-нибудь это пригодится. Во всяком случае в нашем арсенале появится еще один инструмент.

А начнем мы с того, что кратко рассмотрим, как можно получить точки перечения двух кривых другими методами.

Получение точек пересечения кривых аппроксимацией

pascal
   pl := p0; 
   for i := 1 to 20 do begin
      t := i/m;
      pn := PointOnCurve(p0, p1, p2, p3, t);

      pt1 := IntersectBezierLine(
                 pa, pb, pc, pd, 
                 pl, pn);


      pl := pn;
   end;

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

Схема такого решения приведена слева. Здесь p0, p1, p2, p3 - опорные точки первой кривой, pa, pb, pc, pd - опорные точки второй кривой. Мы разбиваем кривую на 20 отрезков. Подставляя для концов отрезков параметр t (который, как известно меняется от нуля до единицы) в уравнение кривой (вызов PointOnCurve), получаем очередную точку на кривой. Для очередного отрезка (pl, pn) вычисляем точку его пересечения со второй кривой (вызов IntersectBezierLine) - pt1. Считаем, что отрезки достаточно малы и пересекают вторую кривую не более чем в одной точке. Отрезки должны быть малы еще и для того, чтобы точнее аппроксимировать первую кривую. То есть чем больше точек мы возьмем, тем точнее получим точки пересечения, но, с другой стороны, придется сделать больше расчетов, что, естественно, скажется на скорости выполнения. В идеале нужно делить первую кривую на отрезки разной длины в зависимости от ее кривизны. Отсюда второй вариант аппроксимации кривой отрезками.

-
-
Кривая Безье очень просто делится пополам. Посмотрите на рисунок слева. Нужно всего лишь поделить отрезки между опорными точками несколько раз пополам и получим две кривые Безье и новые опорные точки. Каждую из кривых можно снова поделить пополам и так далее. После нескольких делений оставшаяся часть кривой будет достаточно близка к отрезку. Собственно из-за этого свойства кривых мы их и используем. Как же определить когда кривая будет достаточно близка к отрезку ? Для наших целей, например, можно найти отклонение средней точки кривой от прямой, проходящей через ее начало и конц. Когда это расстояние станет достаточно малым, деление можно прекратить и искать точки пересечения отрезка со второй кривой.


Таким образом, у нас есть даже два варианта для простого поиска точек пересечения двух кривых Безье. В большинстве случаев они прекрасно работают и дают достаточно близкое приближение, но попробуйте найти точки пересечения, когда кривые почти сливаются, как на рисунке слева.

Теперь перейдем к аналитическому решению.

Вывод системы уравнений

Как известно, кубическая кривая Безье описывается следующим уравнением:

  b(t) = (1-t)3p0 + 3t(1-t)2p1 + 3t2(1-t)p2 + t3p3 

Раскрыв скобки получим полином

  b(t) = At3 + Bt2 + Ct + D 

где

  A = -p0 + 3p1 - 3p2 + p3 
  B = 3p0 - 6p1 + 3p2 
  C = -3p0 + 3p1 
  D = p0 

Таким образом для двух кривых имеем два полинома и параметры для точек пересечения - это корни t1 и t2 системы уравнений

  Axt13 + Bxt12 + Cxt1 + Dx = Ext23 + Fxt22 + Gxt2 + Hx 
  Ayt13 + Byt12 + Cyt1 + Dy = Eyt23 + Fyt22 + Gyt2 + Hy 

Перепишем систему уравнений немного по-другому, заменив переменные на x и y и перенеся все члены в левую часть

  A1x3 + A2x2 + A3x + A4y3 + A5y2 + A6y + A7 = 0 
  B1x3 + B2x2 + B3x + B4y3 + B5y2 + B6y + B7 = 0 

Вот эту систему нам и предстоит решить. Ключ к решению этой системы дает нам теория вычетов[2].

comment
Еще один способ - алгоритм Бухбергера[3] получения редуцированного базиса Гребнера.

Применение теории вычетов для решения системы уравнений

Рассмотрим исходные уравнения системы как полиномы по переменной x .

  f(x) = A1x3 + A2x2 + A3x + A0 
  g(x) = B1x3 + B2x2 + B3x + B0 

здесь

  A0 = A4y3 + A5y2 + A6y + A7 
  B0 = B4y3 + B5y2 + B6y + B7 

Решение системы полиномов от одной переменной состоит в вычислении результанта полиномов. В нашем случае это полиномы f(x) и g(x) . Составим для полиномов f(x) и g(x) матрицу Сильвестра :

$$ S(f,g) = \left( 
   \begin 
       A_1 & A_2 & A_3 & A_0 & 0   & 0   \\ 
       0   & A_1 & A_2 & A_3 & A_0 & 0   \\ 
       0   & 0   & A_1 & A_2 & A_3 & A_0 \\ 
       B_1 & B_2 & B_3 & B_0 & 0   & 0   \\ 
       0   & B_1 & B_2 & B_3 & B_0 & 0   \\ 
       0   & 0   & B_1 & B_2 & B_3 & B_0 \\ 
   \end 
   \right) 
$$

Определитель этой матрицы будет результантом полиномов f(x) и g(x) . Выразим его в символьном виде[4]*:

 det(S(f,g)) = 
 a1^3*b4^3           + a1^3*b0^2*b2        - 2*a1^3*b0*b3*b4     - a1^2*a2*b3*b4^2     - 
 a1^2*a2*b0^2*b1     + a1^2*a2*b0*b3^2     + a1^2*a2*b0*b2*b4    + a1^2*a3*b3^2*b4     - 
 2*a1^2*a3*b2*b4^2   + 2*a1^2*a3*b0*b1*b4  - a1^2*a4*b3^3        + 3*a1^2*a4*b2*b3*b4  - 
 3*a1^2*a4*b1*b4^2   - 2*a1^2*a4*b0*b2^2   + 2*a1^2*a4*b0*b1*b3  + a1*a2^2*b2*b4^2     - 
 a1*a2^2*b0*b2*b3    - a1*a2^2*b0*b1*b4    - a1*a2*a3*b2*b3*b4   + a1*a2*a3*b0*b2^2    + 
 3*a1*a2*a3*b1*b4^2  - 2*a1*a2*a4*b2^2*b4  + a1*a2*a4*b2*b3^2    + 3*a1*a2*a4*b0*b1*b2 + 
 a1*a3^2*b2^2*b4     - 2*a1*a3^2*b1*b3*b4  - a1*a3*a4*b2^2*b3    + 2*a1*a3*a4*b1*b3^2  - 
 2*a1*a3*a4*b0*b1^2  + a1*a4^2*b2^3        + 3*a1*a4^2*b1^2*b4   - 3*a1*a4^2*b1*b2*b3  + 
 a0^2*a1*b1^2*b2     - a0^2*a2*b1^3        + a0*a1^2*b2^2*b4     - a0*a1^2*b2*b3^2     + 
 2*a0*a1^2*b1*b3*b4  - 2*a0*a1^2*b0*b1*b2  + a0*a1*a2*b2^2*b3    - 3*a0*a1*a2*b1*b2*b4 + 
 2*a0*a1*a2*b0*b1^2  - a0*a1*a3*b2^3       - 2*a0*a1*a3*b1^2*b4  + 2*a0*a1*a3*b1*b2*b3 - 
 a2^3*b1*b4^2        + a2^3*b0*b1*b3       + a2^2*a3*b1*b3*b4    - a2^2*a3*b0*b1*b2    - 
 a2^2*a4*b1*b3^2     + 2*a2^2*a4*b1*b2*b4  - 2*a1*a2*a3*b0*b1*b3 - a2^2*a4*b0*b1^2     - 
 a2*a3^2*b1*b2*b4    + a2*a3^2*b0*b1^2     - 3*a2*a3*a4*b1^2*b4  + a2*a3*a4*b1*b2*b3   - 
 a1*a2*a4*b1*b3*b4   + 2*a2*a4^2*b1^2*b3   - a2*a4^2*b1*b2^2     + 2*a0*a2^2*b1^2*b4   - 
 a0*a2^2*b1*b2*b3    + a0*a2*a3*b1*b2^2    + a3^3*b1^2*b4        - a3^2*a4*b1^2*b3     + 
 a3*a4^2*b1^2*b2     - a0*a3^2*b1^2*b2     + a1*a3*a4*b1*b2*b4   - a4^3*b1^3           - 
 2*a0*a1*a4*b1^2*b3  + 2*a0*a3*a4*b1^3     + a0*a1*a4*b1*b2^2    - a0*a2*a4*b1^2*b2 
*Здесь символ ^ означает возведение в степень.

Теперь подставим в это выражение значения A0 и B0 и получим жутковатого вида полином девятой степени:

 h(y) =                                                                                  
                                                                                         
 y^9*(  a1^3*b4^3               - a4^3*b1^3           - 3*a1^2*a4*b1*b4^2                
      + 3*a1*a4^2*b1^2*b4                                                   ) +          
 y^8*(  3*a1^3*b4^2*b5          - 3*a4^2*a5*b1^3      - 6*a1^2*a4*b1*b4*b5               
      - 3*a1^2*a5*b1*b4^2       + 3*a1*a4^2*b1^2*b5   + 6*a1*a4*a5*b1^2*b4  ) +          
 y^7*(  3*a1^3*b4^2*b6          + 3*a1^3*b4*b5^2      - 3*a4^2*a6*b1^3                   
      - 3*a4*a5^2*b1^3          - 6*a1^2*a4*b1*b4*b6  - 3*a1^2*a4*b1*b5^2                
      - 6*a1^2*a5*b1*b4*b5      - 3*a1^2*a6*b1*b4^2   + 3*a1*a4^2*b1^2*b6                
      + 6*a1*a4*a5*b1^2*b5      + 6*a1*a4*a6*b1^2*b4  + 3*a1*a5^2*b1^2*b4   ) +          
 y^6*(  3*a1^3*b4^2*b7          + 6*a1^3*b4*b5*b6     + a1^3*b5^3                        
      - a1^2*a2*b3*b4^2         - 2*a1^2*a3*b2*b4^2   - 3*a4^2*a7*b1^3                   
      - 6*a4*a5*a6*b1^3         - a5^3*b1^3           - 6*a1^2*a4*b1*b4*b7               
      - 6*a1^2*a4*b1*b5*b6      + 3*a1^2*a4*b2*b3*b4  - 6*a1^2*a5*b1*b4*b6               
      - 3*a1^2*a5*b1*b5^2       - 6*a1^2*a6*b1*b4*b5  - 3*a1^2*a7*b1*b4^2                
      + a1*a4^2*b2^3            + 3*a1*a4^2*b1^2*b7   - 3*a1*a4^2*b1*b2*b3               
      - a2^3*b1*b4^2            + 2*a2^2*a4*b1*b2*b4  + 2*a2*a4^2*b1^2*b3                
      - a2*a4^2*b1*b2^2         + a3*a4^2*b1^2*b2     + 6*a1*a4*a5*b1^2*b6               
      + 6*a1*a4*a6*b1^2*b5      + 3*a1*a5^2*b1^2*b5   + 6*a1*a4*a7*b1^2*b4               
      + 6*a1*a5*a6*b1^2*b4      + a1*a2^2*b2*b4^2     - 2*a1*a2*a4*b2^2*b4               
      - 3*a2*a3*a4*b1^2*b4      - a1*a2*a4*b1*b3*b4   + a1*a3*a4*b1*b2*b4                
      + 3*a1*a2*a3*b1*b4^2                                                  ) +          
 y^5*(  6*a1^3*b4*b5*b7         + 3*a1^3*b4*b6^2      + 3*a1^3*b5^2*b6                   
      - 2*a1^2*a2*b3*b4*b5      - 4*a1^2*a3*b2*b4*b5  - 6*a4*a5*a7*b1^3                  
      - 3*a4*a6^2*b1^3          - 3*a5^2*a6*b1^3      - 6*a1^2*a4*b1*b5*b7               
      - 3*a1^2*a4*b1*b6^2       + 3*a1^2*a4*b2*b3*b5  - 6*a1^2*a5*b1*b4*b7               
      - 6*a1^2*a5*b1*b5*b6      + 3*a1^2*a5*b2*b3*b4  - 6*a1^2*a6*b1*b4*b6               
      - 3*a1^2*a6*b1*b5^2       - 6*a1^2*a7*b1*b4*b5  + 2*a1*a4*a5*b2^3                  
      + 6*a1*a4*a5*b1^2*b7      - 6*a1*a4*a5*b1*b2*b3 - 2*a2^3*b1*b4*b5                  
      + 2*a2^2*a4*b1*b2*b5      + 2*a2^2*a5*b1*b2*b4  + 4*a2*a4*a5*b1^2*b3               
      - 2*a2*a4*a5*b1*b2^2      + 2*a3*a4*a5*b1^2*b2  + 6*a1*a4*a6*b1^2*b6               
      + 3*a1*a5^2*b1^2*b6       + 6*a1*a4*a7*b1^2*b5  + 6*a1*a5*a6*b1^2*b5               
      + 2*a1*a2^2*b2*b4*b5      - 2*a1*a2*a4*b2^2*b5  - 3*a2*a3*a4*b1^2*b5               
      -   a1*a2*a4*b1*b3*b5     +   a1*a3*a4*b1*b2*b5 + 6*a1*a5*a7*b1^2*b4               
      + 3*a1*a6^2*b1^2*b4       - 2*a1*a2*a5*b2^2*b4  - 3*a2*a3*a5*b1^2*b4               
      -   a1*a2*a5*b1*b3*b4     +   a1*a3*a5*b1*b2*b4 + 6*a1*a2*a3*b1*b4*b5 )+           
 y^4(   6*a1^3*b4*b6*b7         + 3*a1^3*b5^2*b7      + 3*a1^3*b5*b6^2                   
      - 2*a1^2*a2*b3*b4*b6      -   a1^2*a2*b3*b5^2   - 4*a1^2*a3*b2*b4*b6               
      - 2*a1^2*a3*b2*b5^2       - 6*a4*a6*a7*b1^3     - 3*a5^2*a7*b1^3                   
      - 3*a5*a6^2*b1^3          - 6*a1^2*a4*b1*b6*b7  + 3*a1^2*a4*b2*b3*b6               
      - 6*a1^2*a5*b1*b5*b7      - 3*a1^2*a5*b1*b6^2   + 3*a1^2*a5*b2*b3*b5               
      - 6*a1^2*a6*b1*b4*b7      - 6*a1^2*a6*b1*b5*b6  + 3*a1^2*a6*b2*b3*b4               
      - 6*a1^2*a7*b1*b4*b6      - 3*a1^2*a7*b1*b5^2   + 2*a1*a4*a6*b2^3                  
      +   a1*a5^2*b2^3          + 6*a1*a4*a6*b1^2*b7  - 6*a1*a4*a6*b1*b2*b3              
      + 3*a1*a5^2*b1^2*b7       - 3*a1*a5^2*b1*b2*b3  - 2*a2^3*b1*b4*b6                  
      -   a2^3*b1*b5^2          + 2*a2^2*a4*b1*b2*b6  + 2*a2^2*a5*b1*b2*b5               
      + 2*a2^2*a6*b1*b2*b4      + 4*a2*a4*a6*b1^2*b3  + 2*a2*a5^2*b1^2*b3                
      - 2*a2*a4*a6*b1*b2^2      -   a2*a5^2*b1*b2^2   + 2*a3*a4*a6*b1^2*b2               
      +   a3*a5^2*b1^2*b2       + 6*a1*a4*a7*b1^2*b6  + 6*a1*a5*a6*b1^2*b6               
      + 2*a1*a2^2*b2*b4*b6      +   a1*a2^2*b2*b5^2   - 2*a1*a2*a4*b2^2*b6               
      - 3*a2*a3*a4*b1^2*b6      -   a1*a2*a4*b1*b3*b6 +   a1*a3*a4*b1*b2*b6              
      + 6*a1*a5*a7*b1^2*b5      + 3*a1*a6^2*b1^2*b5   - 2*a1*a2*a5*b2^2*b5               
      - 3*a2*a3*a5*b1^2*b5      -   a1*a2*a5*b1*b3*b5 +   a1*a3*a5*b1*b2*b5              
      + 6*a1*a6*a7*b1^2*b4      + 6*a1*a2*a3*b1*b4*b6 + 3*a1*a2*a3*b1*b5^2               
      - 2*a1*a2*a6*b2^2*b4      - 3*a2*a3*a6*b1^2*b4  -   a1*a2*a6*b1*b3*b4              
      +   a1*a3*a6*b1*b2*b4                                                 )+           
 y^3*(  3*a1^3*b4*b7^2          + 6*a1^3*b5*b6*b7     +   a1^3*b6^3                      
      - 2*a1^2*a2*b3*b4*b7      - 2*a1^2*a2*b3*b5*b6  - 4*a1^2*a3*b2*b4*b7               
      - 4*a1^2*a3*b2*b5*b6      +   a1^2*a3*b3^2*b4   - 3*a4*a7^2*b1^3                   
      - 6*a5*a6*a7*b1^3         -   a6^3*b1^3         -   a1^2*a4*b3^3                   
      - 3*a1^2*a4*b1*b7^2       + 3*a1^2*a4*b2*b3*b7  - 6*a1^2*a5*b1*b6*b7               
      + 3*a1^2*a5*b2*b3*b6      - 6*a1^2*a6*b1*b5*b7  - 3*a1^2*a6*b1*b6^2                
      + 3*a1^2*a6*b2*b3*b5      - 6*a1^2*a7*b1*b4*b7  - 6*a1^2*a7*b1*b5*b6               
      + 3*a1^2*a7*b2*b3*b4      + 2*a1*a4*a7*b2^3     + 2*a1*a5*a6*b2^3                  
      + 6*a1*a4*a7*b1^2*b7      - 6*a1*a4*a7*b1*b2*b3 + 6*a1*a5*a6*b1^2*b7               
      - 6*a1*a5*a6*b1*b2*b3     - 2*a2^3*b1*b4*b7     - 2*a2^3*b1*b5*b6                  
      -   a2^2*a4*b1*b3^2       + 2*a2^2*a4*b1*b2*b7  + 2*a2^2*a5*b1*b2*b6               
      + 2*a2^2*a6*b1*b2*b5      +   a2^2*a3*b1*b3*b4  + 2*a2^2*a7*b1*b2*b4               
      + 4*a2*a4*a7*b1^2*b3      + 4*a2*a5*a6*b1^2*b3  - 2*a2*a4*a7*b1*b2^2               
      - 2*a2*a5*a6*b1*b2^2      +   a3^3*b1^2*b4      -   a3^2*a4*b1^2*b3                
      + 2*a3*a4*a7*b1^2*b2      + 2*a3*a5*a6*b1^2*b2  + 2*a1*a2^2*b2*b4*b7               
      + 2*a1*a2^2*b2*b5*b6      - 2*a1*a2*a4*b2^2*b7  +   a1*a2*a4*b2*b3^2               
      +   a1*a3^2*b2^2*b4       - 2*a1*a3^2*b1*b3*b4  -   a1*a3*a4*b2^2*b3               
      + 2*a1*a3*a4*b1*b3^2      -   a2*a3^2*b1*b2*b4  - 3*a2*a3*a4*b1^2*b7               
      +   a2*a3*a4*b1*b2*b3     -   a1*a2*a4*b1*b3*b7 +   a1*a3*a4*b1*b2*b7              
      + 6*a1*a5*a7*b1^2*b6      + 3*a1*a6^2*b1^2*b6   - 2*a1*a2*a5*b2^2*b6               
      - 3*a2*a3*a5*b1^2*b6      -   a1*a2*a5*b1*b3*b6 +   a1*a3*a5*b1*b2*b6              
      + 6*a1*a6*a7*b1^2*b5      + 3*a1*a7^2*b1^2*b4   + 6*a1*a2*a3*b1*b4*b7              
      + 6*a1*a2*a3*b1*b5*b6     -   a1*a2*a3*b2*b3*b4 - 2*a1*a2*a6*b2^2*b5               
      - 2*a1*a2*a7*b2^2*b4      - 3*a2*a3*a6*b1^2*b5  - 3*a2*a3*a7*b1^2*b4               
      -   a1*a2*a6*b1*b3*b5     -   a1*a2*a7*b1*b3*b4 +   a1*a3*a6*b1*b2*b5              
      +   a1*a3*a7*b1*b2*b4                                                 )+           
 y^2*(+ 3*a1^3*b5*b7^2          + 3*a1^3*b6^2*b7      - 2*a1^2*a2*b3*b5*b7               
      -   a1^2*a2*b3*b6^2       - 4*a1^2*a3*b2*b5*b7  - 2*a1^2*a3*b2*b6^2                
      +   a1^2*a3*b3^2*b5       - 3*a5*a7^2*b1^3      - 3*a6^2*a7*b1^3                   
      -   a1^2*a5*b3^3          - 3*a1^2*a5*b1*b7^2   + 3*a1^2*a5*b2*b3*b7               
      - 6*a1^2*a6*b1*b6*b7      + 3*a1^2*a6*b2*b3*b6  - 6*a1^2*a7*b1*b5*b7               
      - 3*a1^2*a7*b1*b6^2       + 3*a1^2*a7*b2*b3*b5  + 2*a1*a5*a7*b2^3                  
      +   a1*a6^2*b2^3          + 6*a1*a5*a7*b1^2*b7  - 6*a1*a5*a7*b1*b2*b3              
      + 3*a1*a6^2*b1^2*b7       - 3*a1*a6^2*b1*b2*b3  - 2*a2^3*b1*b5*b7                  
      -   a2^3*b1*b6^2          -   a2^2*a5*b1*b3^2   + 2*a2^2*a5*b1*b2*b7               
      + 2*a2^2*a6*b1*b2*b6      +   a2^2*a3*b1*b3*b5  + 2*a2^2*a7*b1*b2*b5               
      + 4*a2*a5*a7*b1^2*b3      + 2*a2*a6^2*b1^2*b3   - 2*a2*a5*a7*b1*b2^2               
      -   a2*a6^2*b1*b2^2       +   a3^3*b1^2*b5      -   a3^2*a5*b1^2*b3                
      + 2*a3*a5*a7*b1^2*b2      +   a3*a6^2*b1^2*b2   + 2*a1*a2^2*b2*b5*b7               
      +   a1*a2^2*b2*b6^2       - 2*a1*a2*a5*b2^2*b7  +   a1*a2*a5*b2*b3^2               
      +   a1*a3^2*b2^2*b5       - 2*a1*a3^2*b1*b3*b5  -   a1*a3*a5*b2^2*b3               
      + 2*a1*a3*a5*b1*b3^2      -   a2*a3^2*b1*b2*b5  - 3*a2*a3*a5*b1^2*b7               
      +   a2*a3*a5*b1*b2*b3     -   a1*a2*a5*b1*b3*b7 +   a1*a3*a5*b1*b2*b7              
      + 6*a1*a6*a7*b1^2*b6      + 3*a1*a7^2*b1^2*b5   + 6*a1*a2*a3*b1*b5*b7              
      + 3*a1*a2*a3*b1*b6^2      -   a1*a2*a3*b2*b3*b5 - 2*a1*a2*a6*b2^2*b6               
      - 2*a1*a2*a7*b2^2*b5      - 3*a2*a3*a6*b1^2*b6  - 3*a2*a3*a7*b1^2*b5               
      -   a1*a2*a6*b1*b3*b6     -   a1*a2*a7*b1*b3*b5 +   a1*a3*a6*b1*b2*b6              
      +   a1*a3*a7*b1*b2*b5                                                 )+           
 y*(  + 3*a1^3*b6*b7^2          - 2*a1^2*a2*b3*b6*b7  - 4*a1^2*a3*b2*b6*b7               
      +   a1^2*a3*b3^2*b6       - 3*a6*a7^2*b1^3      -   a1^2*a6*b3^3                   
      - 3*a1^2*a6*b1*b7^2       + 3*a1^2*a6*b2*b3*b7  - 6*a1^2*a7*b1*b6*b7               
      + 3*a1^2*a7*b2*b3*b6      + 2*a1*a6*a7*b2^3     + 6*a1*a6*a7*b1^2*b7               
      - 6*a1*a6*a7*b1*b2*b3     + 3*a1*a7^2*b1^2*b6   - 2*a2^3*b1*b6*b7                  
      -   a2^2*a6*b1*b3^2       + 2*a2^2*a6*b1*b2*b7  +   a2^2*a3*b1*b3*b6               
      + 2*a2^2*a7*b1*b2*b6      + 4*a2*a6*a7*b1^2*b3  - 2*a2*a6*a7*b1*b2^2               
      +   a3^3*b1^2*b6          -   a3^2*a6*b1^2*b3   + 2*a3*a6*a7*b1^2*b2               
      + 2*a1*a2^2*b2*b6*b7      + 6*a1*a2*a3*b1*b6*b7 -   a1*a2*a3*b2*b3*b6              
      - 2*a1*a2*a6*b2^2*b7      +   a1*a2*a6*b2*b3^2  - 2*a1*a2*a7*b2^2*b6               
      +   a1*a3^2*b2^2*b6       - 2*a1*a3^2*b1*b3*b6  -   a1*a3*a6*b2^2*b3               
      + 2*a1*a3*a6*b1*b3^2      -   a2*a3^2*b1*b2*b6  - 3*a2*a3*a6*b1^2*b7               
      +   a2*a3*a6*b1*b2*b3     - 3*a2*a3*a7*b1^2*b6  -   a1*a2*a6*b1*b3*b7              
      -   a1*a2*a7*b1*b3*b6     +   a1*a3*a6*b1*b2*b7 +   a1*a3*a7*b1*b2*b6 )+           
 (                                                                                       
      + a1^3*b7^3               - a1^2*a2*b3*b7^2     - 2*a1^2*a3*b2*b7^2                
      + a1^2*a3*b3^2*b7         - a7^3*b1^3           - a1^2*a7*b3^3                     
      - 3*a1^2*a7*b1*b7^2       + 3*a1^2*a7*b2*b3*b7  + a1*a7^2*b2^3                     
      + 3*a1*a7^2*b1^2*b7       - 3*a1*a7^2*b1*b2*b3  - a2^3*b1*b7^2                     
      + a2^2*a3*b1*b3*b7        - a2^2*a7*b1*b3^2     + 2*a2^2*a7*b1*b2*b7               
      + 2*a2*a7^2*b1^2*b3       - a2*a7^2*b1*b2^2     + a3^3*b1^2*b7                     
      - a3^2*a7*b1^2*b3         + a3*a7^2*b1^2*b2     + a1*a2^2*b2*b7^2                  
      + 3*a1*a2*a3*b1*b7^2      - a1*a2*a3*b2*b3*b7   - 2*a1*a2*a7*b2^2*b7               
      + a1*a2*a7*b2*b3^2        + a1*a3^2*b2^2*b7     - 2*a1*a3^2*b1*b3*b7               
      - a1*a3*a7*b2^2*b3        + 2*a1*a3*a7*b1*b3^2  - a2*a3^2*b1*b2*b7                 
      - 3*a2*a3*a7*b1^2*b7      + a2*a3*a7*b1*b2*b3   - a1*a2*a7*b1*b3*b7                
      + a1*a3*a7*b1*b2*b7                                                   )            

Мы можем найти корни уравнения h(y) = 0 , например, используя метод Ньютона (метод касательных. Здесь мы не будем на нем останавливаться). Полученные корни уравнения y подставим в одно из уравнений системы и, решив кубическое уравнение, получим значения x .

Вырожденные случаи

-
-
Обратили внимание, что мы ничего не сказали о случаях, когда результант равен нулю ? Например, если коэффициенты A1 и B1 одновременно будут равны нулю, то и определитель матрицы и все коэффициенты нашего полинома также будут равны нулю.

Например, кривые Безье, заданные точками:

 (0, 0) (0, 2) (1, 3) (3, 3)
 (0, 4) (0, 2) (1, 1) (3, 1)

описываются следующими уравнениями:

 3x2 - 3y2 = 0
 -3x2 + 6x - 3y2 + 6y - 4 = 0

Для такой системы уравнение h(y) =0 вырождается в равенство 0 = 0 . Между тем, система имеет очевидное решение ( x = y ) и его легко получить, решив квадратное уравнение.

Давайте выделим такой случай в отдельное решение. Начнем сначала. Полиномы f(x) и g(x) имеют вид:

  f(x) = A2x2 + A3x + A0 
  g(x) = B2x2 + B3x + B0 

Составим матрицу Сильвестра для этих полиномов

$$ S(f,g) = \left( 
   \begin 
       A_2 & A_3 & A_0 & 0     \\ 
       0   & A_2 & A_3 & A_0   \\ 
       B_2 & B_3 & B_0 & 0     \\ 
       0   & B_2 & B_3 & B_0   \\ 
   \end 
   \right) 
$$

Найдем определитель этой системы:

  det(S(f,g)) =                   
    a2^2*b0^2      - a2*a3*b0*b3  + a0^2*b2^2      + a0*a2*b3^2   
  - 2*a0*a2*b0*b2  + a3^2*b0*b2   - a0*a3*b2*b3                   

И, наконец, подставив значения A0 и B0 в это выражение, получим полином шестой степени:

  h(y) =   
                                                                                 
  y^6*( a2^2*b4^2     + a4^2*b2^2     - 2*a2*a4*b2*b4                      ) +   
  y^5*( 2*a2^2*b4*b5  + 2*a4*a5*b2^2  - 2*a2*a4*b2*b5  - 2*a2*a5*b2*b4     ) +   
  y^4*( 2*a2^2*b4*b6  + a2^2*b5^2     + 2*a4*a6*b2^2   + a5^2*b2^2     -         
        2*a2*a4*b2*b6 - 2*a2*a5*b2*b5 - 2*a2*a6*b2*b4                      ) +    
  y^3*( 2*a2^2*b4*b7  + 2*a2^2*b5*b6  - a2*a3*b3*b4    + 2*a4*a7*b2^2  +         
        2*a5*a6*b2^2  + a2*a4*b3^2    - 2*a2*a4*b2*b7  + a3^2*b2*b4    -         
        a3*a4*b2*b3   - 2*a2*a5*b2*b6 - 2*a2*a6*b2*b5  - 2*a2*a7*b2*b4     ) +   
  y^2*( 2*a2^2*b5*b7  + a2^2*b6^2     - a2*a3*b3*b5    + 2*a5*a7*b2^2  +         
        a6^2*b2^2     + a2*a5*b3^2    - 2*a2*a5*b2*b7  + a3^2*b2*b5    -         
        a3*a5*b2*b3   - 2*a2*a6*b2*b6 - 2*a2*a7*b2*b5                      ) +   
  y*(   2*a2^2*b6*b7  - a2*a3*b3*b6   + 2*a6*a7*b2^2   + a2*a6*b3^2    -         
        2*a2*a6*b2*b7 - 2*a2*a7*b2*b6 + a3^2*b2*b6     - a3*a6*b2*b3       ) +   
    (   a2^2*b7^2     - a2*a3*b3*b7   + a7^2*b2^2      + a2*a7*b3^2    -         
        2*a2*a7*b2*b7 + a3^2*b2*b7    - a3*a7*b2*b3                        )     

Точка самопересечения кривой Безье

-
-
Чтобы получить точку самопересечения кривой, нужно решить систему:

  A1t13 + A2t12 + A3t1 = A1t23 + A2t22 + A3t2 
  B1t13 + B2t12 + B3t1 = B1t23 + B2t22 + B3t2, t1 ≠ t2 

Примем:

  t2 = t1 + a, a ≠ 0 

Подставив это выражение в исходную систему и сократив на a , получим следующую систему уравнений:

  3A1t12 + (3A1a + 2A2)t1 + (A1a2 + A2a + A3) = 0 
  3B1t12 + (3B1a + 2B2)t1 + (B1a2 + B2a + B3) = 0 

Рассматривая полиномы системы по переменной t получим результант. Подставив в него коэффициенты, получим следующее уравнение для переменной a :

 C = (3*a1^2*b2^2     - 6*a1*a2*b1*b2  + 3*a2^2*b1^2)
 D = (9*a1^2*b3^2     - 12*a1*a2*b2*b3 + 12*a1*a3*b2^2  -
      18*a1*a3*b1*b3  + 12*a2^2*b1*b3  - 12*a2*a3*b1*b2 + 9*a3^2*b1^2)
 a^2 = -D/C 

Если это уравнение имеет решение, то его можно подставить в первое уравнение системы и, решив квадратное уравнение, получить значение t1 .

Вместо заключения

-
-
На этом, пожалуй, остановимся. Мы получили аналитическое решение и вкратце познакомились с теорией вычетов.

1Kevin Lindsey JavaScript, посетите также сайт KevLinDev
3Д.Кокс, Дж.Литтл, Д.О'Ши Идеалы, многообразия и алгоритмы.
4Ни один из математических пакетов не пострадал не использовался для получения выражений в этой статье.

04 ноября 2010

Comments
27.11.2010 03:02
 На псевдокоде это бы воспринималось бы лучше нежели голая математика :)
 Кстати "колесо ссылок" на мозиле не пашет.
 
27.11.2010 11:23 alexBlack
 По-моему здесь псевдокод ни к чему. Решение получено, а как пример реализации 
 можно посмотреть java-скрипт по первой ссылке.

 С "колесом" посмотрим. Спасибо за замечание. 
 
21.12.2010 03:51 iMerlin
 Вэри тханькс!!! :-)
 Буду изучать - как раз ко времени.
 
29.04.2011 11:58 seoblogst
 avtor molodec!
04.05.2011 02:00 seomaestro
 Шикарный пост! Очень понравился)
 
04.05.2011 08:43 vevlohenor
 Пост что надо..
 
10.05.2011 05:00 seonickname
 Хороший пост, пожалуй ретвитну)))
 
11.05.2011 01:52 seoblya
 100% подпишусь на Вас по RSS

20.05.2011 08:49 unmaspJamyday
 Hi !!! Good job! 
 
28.05.2011 08:50 eve isk
      
 I've been looking for this info for a long time, thanks for your work.

05.06.2011 10:52 Newautomen
 Скажите, у Вас есть рассылка новостей на сайбскрайбе или маиллисте, 
 просто интернет ограничен, было бы удобно получать материал на почту.
 
06.06.2011 08:45 alexBlack
 Пока ничего подобного нет, да и новостей немного. 
 Если нужно, напишите мне на e-mail, буду сообщать по мере обновлений.
 (e-mail есть на страничке контактов) 
 
12.06.2011 02:33 slerryDem
 I think this post is concerning to study. Keep it going guy
 
19.06.2011 12:45 bayan
 Спасибо, полезный материал. Добавил в закладки.

12.07.2016 05:14 pro_100_gram
 Спасибо большое! Очень полезный материал.
 
Вы можете оставить комментарий или задать вопрос
Ваше имя:

Текст сообщения:


Copyright © 2009-2014 by