Если требуется превратить корявую ломаную в симпатичную гладкую кривую так, чтобы она проходила через все точки ломаной, то тут может помочь кубический сплайн Эрмита. Если лень читать теорию, а её здесь много, то сразу идём в конец статьи, там исходники и интерактив.
Интерполяция сплайнами простыми словами
Интерполяция сплайнами – это умный способ провести идеально гладкую линию через набор точек. Вместо одной корявой кривой на всех точках или ломаной, она использует много коротких простых кривых, которые магически склеиваются в точках так, что линия получается плавной, без углов и изломов, и при этом точно проходит через каждую заданную точку.
Что такое сплайн
Идея сплайна — вместо одной сложной кривой используем много простых. Представим себе гибкую деревянную или пластиковую линейку (именно такую линейку раньше называли «сплайном», от английского «spline» – рейка, планка). Раньше чертёжники брали такую гибкую рейку, пришпиливали ее гвоздиками в нужных местах, и рейка сама изгибалась, образуя плавную, минимально напряженную кривую между гвоздиками. Физика заставляла ее быть гладкой и естественной.
Математический сплайн – это цифровой аналог этой гибкой рейки. Мы разбиваем наш набор точек на отрезки (между соседними точками — гвоздиками). На каждом таком маленьком отрезке мы строим свою простую плавную кривую (чаще всего – параболу или кубическую кривую). Эти маленькие кривые «склеиваются» друг с другом в наших точках (гвоздиках) очень аккуратно.
Склейка
Самое важное в сплайнах – как эти маленькие кривые соединяются в точках. Мы требуем не просто, чтобы они пересекались в точке (это очевидно). Мы требуем, чтобы в точке соединения:
- Не было «угла»: Кривые подходят к точке под одинаковым углом (совпадают первые производные). Это гарантирует гладкость, как будто это одна кривая.
- Не было «излома»: Кривизна линии меняется плавно при переходе через точку (совпадают вторые производные). Это убирает резкие перегибы.
Почему кубические
Чаще всего используют кубические сплайны (строят на каждом отрезке между точками маленькую кубическую кривую):
|
1 |
y = a + b*x + c*x^2 + d*x^3 |
Кубики – это «золотая середина»:
- Они достаточно простые для вычислений.
- Они обладают достаточной гибкостью, чтобы изображать разные изгибы (прямые, S-образные кривые).
- При их склейке мы можем обеспечить совпадение и значения, и наклона, и кривизны (как раз те самые первые и вторые производные), что дает очень хорошую гладкость.
Кубические сплайны
У нас есть n точек данных: (x₀, y₀), (x₁, y₁), ..., (xₙ₋₁, yₙ₋₁). Мы разбиваем область на n-1 интервалов: [x₀, x₁], [x₁, x₂], ..., [xₙ₋₂, xₙ₋₁].
На каждом интервале [xᵢ, xᵢ₊₁] мы строим свою кубическую функцию Sᵢ(x):Sᵢ(x) = aᵢ + bᵢ(x - xᵢ) + cᵢ(x - xᵢ)² + dᵢ(x - xᵢ)³
aᵢ, bᵢ, cᵢ, dᵢ— коэффициенты, которые нам нужно найти для каждого сегментаi.(x - xᵢ)— это смещение внутри текущего интервала. Так запись часто удобнее.
Чтобы вся кривая S(x), составленная из кусочков Sᵢ(x), была гладкой и проходила через точки, мы накладываем четыре типа условий:
Прохождение через точки: Условие Интерполяции
На ЛЕВОМ конце своего интервала [xᵢ, xᵢ₊₁] каждая функция Sᵢ(x) должна проходить через точку (xᵢ, yᵢ):Sᵢ(xᵢ) = yᵢ для всех i = 0, 1, ..., n-2
Подставляя x = xᵢ в уравнение Sᵢ(x):aᵢ + bᵢ(0) + cᵢ(0)² + dᵢ(0)³ = yᵢ => aᵢ = yᵢ
Коэффициент aᵢ — это просто значение yᵢ в левой точке интервала! Магия упрощается.
На ПРАВОМ конце своего интервала [xᵢ, xᵢ₊₁] каждая функция Sᵢ(x) должна проходить через точку (xᵢ₊₁, yᵢ₊₁):Sᵢ(xᵢ₊₁) = yᵢ₊₁ для всех i = 0, 1, ..., n-2
Подставляя x = xᵢ₊₁ (т.е. (x - xᵢ) = hᵢ, где hᵢ = xᵢ₊₁ - xᵢ — длина интервала):aᵢ + bᵢ(hᵢ) + cᵢ(hᵢ)² + dᵢ(hᵢ)³ = yᵢ₊₁
Это дает нам первое уравнение связывающее bᵢ, cᵢ, dᵢ.
Совпадение наклонов в узлах: Условие Непрерывности Первой Производной
Производная Sᵢ'(x) описывает наклон (угол наклона касательной) нашей «рейки» в точке x.
В каждой внутренней точке стыка xᵢ (для i = 1, 2, ..., n-2) наклон, с которым подходит «рейка» слева S'ᵢ₋₁(xᵢ), должен равняться наклону, с которым начинается «рейка» справа S'ᵢ(xᵢ).
Считаем производную кубической функции:Sᵢ'(x) = bᵢ + 2cᵢ(x - xᵢ) + 3dᵢ(x - xᵢ)²
Условие гладкости в точке xᵢ:S'ᵢ₋₁(xᵢ) = S'ᵢ(xᵢ)
Рассматриваем левую функцию (Sᵢ₋₁) в ее правом конце (x = xᵢ, для нее (x - xᵢ₋₁) = hᵢ₋₁):S'ᵢ₋₁(xᵢ) = bᵢ₋₁ + 2cᵢ₋₁(hᵢ₋₁) + 3dᵢ₋₁(hᵢ₋₁)²
Рассматриваем правую функцию (Sᵢ) в ее левом конце (x = xᵢ, для нее (x - xᵢ) = 0):S'ᵢ(xᵢ) = bᵢ + 2cᵢ(0) + 3dᵢ(0)² = bᵢ
Уравнение:bᵢ₋₁ + 2cᵢ₋₁(hᵢ₋₁) + 3dᵢ₋₁(hᵢ₋₁)² = bᵢ для i = 1, 2, ..., n-2
Это дает нам второе уравнение (связывает коэффициенты соседних сегментов).
Совпадение кривизны в узлах: Условие Непрерывности Второй Производной
Вторая производная Sᵢ''(x) описывает скорость изменения наклона, то есть кривизну (насколько резко поворачивает рейка).
В каждой внутренней точке стыка xᵢ (для i = 1, 2, ..., n-2) кривизна, с которой подходит «рейка» слева S''ᵢ₋₁(xᵢ), должна равняться кривизне, с которой начинается «рейка» справа S''ᵢ(xᵢ).
Считаем вторую производную:Sᵢ''(x) = 2cᵢ + 6dᵢ(x - xᵢ)
Условие гладкости в точке xᵢ:S''ᵢ₋₁(xᵢ) = S''ᵢ(xᵢ)
Левая функция (Sᵢ₋₁) в ее правом конце (x = xᵢ, (x - xᵢ₋₁) = hᵢ₋₁):S''ᵢ₋₁(xᵢ) = 2cᵢ₋₁ + 6dᵢ₋₁(hᵢ₋₁)
Правая функция (Sᵢ) в ее левом конце (x = xᵢ, (x - xᵢ) = 0):S''ᵢ(xᵢ) = 2cᵢ + 6dᵢ(0) = 2cᵢ
Уравнение:2cᵢ₋₁ + 6dᵢ₋₁(hᵢ₋₁) = 2cᵢ
Делим обе части на 2:cᵢ₋₁ + 3dᵢ₋₁(hᵢ₋₁) = cᵢ для i = 1, 2, ..., n-2
Это дает нам третье уравнение (тоже связывает коэффициенты соседних сегментов).
Граничные Условия: Что происходит на концах всей «цепочки»?
У нас есть две дополнительные точки (x₀ и xₙ₋₁), где нужно задать поведение. Без этого система не замкнется. Самые распространенные варианты:
Естественный сплайн (Natural Spline):
Кривизна (вторая производная) на концах равна 0. Рейка «выходит» и «входит» прямо.
S''₀(x₀) = 0=>2c₀ + 6d₀(0) = 0=>c₀ = 0S''ₙ₋₂(xₙ₋₁) = 0=>2cₙ₋₂ + 6dₙ₋₂(hₙ₋₂) = 0(т.к. для последнего сегментаx = xₙ₋₁это(x - xₙ₋₂) = hₙ₋₂)
Закрепленный сплайн (Clamped Spline):
Мы задаем конкретные значения наклона (первой производной) на концах. Например, чтобы рейка начиналась и заканчивалась горизонтально.
S'₀(x₀) = заданный_наклон_A=>b₀ + 2c₀(0) + 3d₀(0)² = b₀ = заданный_наклон_AS'ₙ₋₂(xₙ₋₁) = заданный_наклон_B=>bₙ₋₂ + 2cₙ₋₂(hₙ₋₂) + 3dₙ₋₂(hₙ₋₂)² = заданный_наклон_B
Собираем «Пазл»: Система Уравнений:
- Количество неизвестных: У нас
n-1сегментов. На каждый сегмент по 4 коэффициента (aᵢ, bᵢ, cᵢ, dᵢ). Итого4*(n-1)неизвестных. - Количество уравнений:
- Интерполяция (левый конец):
n-1уравнений (по одному на сегмент,aᵢ = yᵢ). - Интерполяция (правый конец):
n-1уравнений (Sᵢ(xᵢ₊₁) = yᵢ₊₁). - Непрерывность 1-й производной:
n-2уравнений (по одному на каждый внутренний узел). - Непрерывность 2-й производной:
n-2уравнений (по одному на каждый внутренний узел). - Граничные условия:
2уравнения (по одному на каждый конец). - Итого:
(n-1) + (n-1) + (n-2) + (n-2) + 2 = 4n - 4уравнений. То есть4*(n-1)уравнений. Ровно столько, сколько неизвестных!
- Интерполяция (левый конец):
- Структура системы: Полученная система линейных алгебраических уравнений (СЛАУ) имеет специальную «ленточную» (или «трехдиагональную») структуру: каждое уравнение связывает только коэффициенты 2-3 соседних сегментов. Это ключевой момент, позволяющий решать ее очень эффективно (методом прогонки или прямыми методами для ленточных матриц), даже для большого числа точек.
Сплайн Эрмита: Суть
Представьте, что вы не просто ставите точки (гвоздики) на доске, но и задаете угол, под которым гибкая рейка должна подходить к каждому гвоздику. То есть вы контролируете не только где проходит кривая, но и как она входит и выходит из каждой точки.
Кубический Сплайн Эрмита: Математическая Магия
Рассмотрим один сегмент между двумя точками: P₀ = (x₀, y₀) и P₁ = (x₁, y₁).
Мы знаем не только положения точек, но и значения первой производной (наклоны) в них:
P'₀ = dy/dxв точкеP₀(обозначимm₀)P'₁ = dy/dxв точкеP₁(обозначимm₁)
Наша цель:
Найти кубический многочлен H(x) = ax³ + bx² + cx + d для интервала [x₀, x₁], удовлетворяющий четырем условиям:
H(x₀) = y₀(Кривая проходит через начальную точку)H(x₁) = y₁(Кривая проходит через конечную точку)H'(x₀) = m₀(Наклон в начальной точке задан)H'(x₁) = m₁(Наклон в конечной точке задан)
Решение: Вывод коэффициентов
Производная:
H'(x) = 3ax² + 2bx + c
Записываем систему уравнений:
- Из условия 1:
a(x₀)³ + b(x₀)² + c(x₀) + d = y₀ - Из условия 2:
a(x₁)³ + b(x₁)² + c(x₁) + d = y₁ - Из условия 3:
3a(x₀)² + 2b(x₀) + c = m₀ - Из условия 4:
3a(x₁)² + 2b(x₁) + c = m₁
Упростим, введя длину сегмента:
h = x₁ - x₀
Решим систему. Решать её громоздко, но результат очень изящен. Коэффициенты можно выразить через известные y₀, y₁, m₀, m₁, h:
a = ( 2(y₀ - y₁) + h(m₀ + m₁) ) / h³b = ( 3(y₁ - y₀) - h(2m₀ + m₁) ) / h²c = m₀d = y₀
Элегантная форма: Базисные Функции Эрмита
Математическую красоту сплайнам Эрмита придает представление через кубические базисные функции Эрмита. Кривая на сегменте выражается как взвешенная сумма этих базисов, умноженных на заданные значения положения и наклона:
H(x) = y₀ * H₀₀(t) + y₁ * H₀₁(t) + m₀ * H₁₀(t) + m₁ * H₁₁(t)
где t = (x - x₀) / h — нормированный параметр (от 0 в x₀ до 1 в x₁).
Сами базисные функции (определенные на [0, 1]):
Вес для начальной позиции y₀: |
H₀₀(t) = (1 - t)²(1 + 2t) = 2t³ - 3t² + 1 |
Вес для конечной позиции y₁: |
H₀₁(t) = t²(3 - 2t) = -2t³ + 3t² |
Вес для начального наклона m₀: |
H₁₀(t) = t(1 - t)² = t³ - 2t² + t |
Вес для конечного наклона m₁: |
H₁₁(t) = t²(t - 1) = t³ - t² |
Сплайн Эрмита: Реализация
Выделим основное из теоретического блока.
Вычисление точек
Мы должны пройтись по точкам массива, взять пару соседних точек Pi и Pi+1 в качестве сегмента, разбить его N-1 частей и выполнить вычисления:
|
1 2 |
pt.X := Hs[j,0]*P1.X + Hs[j,1]*M0.X + Hs[j,2]*P2.X + Hs[j,3]*M1.X; pt.Y := Hs[j,0]*P1.Y + Hs[j,1]*M0.Y + Hs[j,2]*P2.Y + Hs[j,3]*M1.Y; |
Где j-индекс вычисляемой точки внутри сегмента, Hs[j, 0..3] — заранее посчитанные базисные функции Эрмита. P1 и P2 — точки начала и конца рассматриваемого сегмента. M0 и M1 — векторы касательных в начале и конце сегмента.
Базисные функции Эрмита
Базисный функции Эрмита вычисляются только на основании параметра t. Параметр t — это относительное положение рассматриваемого интервала на сегменте. Всегда нормирован от 0 до 1.
Каждый отрезок ломаной делится на одинаковое, заранее заданное, количество частей. Вне зависимости от его длины, при вычислении координат t будет одинаков для каждого интервала внутри рассматриваемого сегмента.
Поэтому, мы кэшируем функции Эрмита заранее, до вычислений.
|
1 2 3 4 5 6 |
t := j / (PointsPerSegment-1); t_t := t*t; Hs[j,0] := t_t*( 2*t - 3) + 1; // H00 = 2t^3 - 3t^2+1 Hs[j,1] := t_t*( t - 2) + t; // H10 = t^3 - 2t^2+t Hs[j,2] := t_t*(-2*t + 3); // H01 =-2t^3 + 3t^2 Hs[j,3] := t_t*( t - 1); // H11 = t^3 - t^2 |
Где j-индекс вычисляемой точки внутри сегмента.
Касательные
Сплайн Эрмита характерен тем, что мы можем задавать входящий и исходящий наклон для каждого сегмента. Возможно, имело бы смысл завернуться на ручную коррекцию наклонов в точках ломаной. Это было бы чрезвычайно интересно, но пока руки не дошли 😞. В реализации мы будем просто учитывать данные точек, находящихся до и после вычисляемого сегмента. Поэтому, для расчётов производных на концах сегмента нам дополнительно понадобятся точки Pi-1 и Pi+2.
Если ломаная незамкнутая, то в качестве предыдущей и следующей точек на концах массива берём крайние точки массива.
Расчёт производных выглядит следующим образом:
Mi = ((Pi+1-Pi) + (Pi-Pi-1))/2 = (Pi+1 - Pi-1)/2
|
1 2 3 |
// Касательные векторы с учетом натяжения M0 := (P2 - P0)*Tension*0.5; M1 := (P3 - P1)*Tension*0.5; |
Где P1 и P2 — точки начала и конца рассматриваемого сегмента. P0 — точка начала предыдущего сегмента. P3 — точка конца следующего сегмента.
Натяжение
Введём такое понятие, как натяжение (Tension). Это некий коэффициент, влияющий на значение наклона. Производные будем домножать на этот коэффициент.
Mi = Tension * (Pi+1 - Pi-1)/2
По уму, надо было бы сделать так, чтобы при значении =0, коэффициент никак не влиял на значения производных. То есть, домножать на самом деле на (1-Tension). Но в таком случае мы будем вынуждены болтаться в интервале [0,1], и это, поверьте, скучно. Ниже, в интерактивной части, можно посмотреть, как этот параметр влияет на форму кривой.
Листинг
|
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 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 |
type TPointFDynArray = array of TPointF; function CreateHermiteSplinesPoints(const Points: TPointFDynArray; Tension: Single = 0.5; PointsPerSegment: Integer = 10; Closed: Boolean = True): TPointFDynArray; var N, L, Idx, i, j: Integer; P0, P1, P2, P3, pt: TPointF; M0, M1: TPointF; t, t_t: Double; Hs: array of array[0..3] of Double; begin N := Length(Points); if (N<2) or (Tension<=0) or (PointsPerSegment<2) then exit(Points); // Формируем кэш базисных функций Эрмита SetLength(Hs, PointsPerSegment); for j := 0 to PointsPerSegment-1 do begin t := j / (PointsPerSegment-1); t_t := t*t; // Базисные функции Эрмита Hs[j,0] := t_t*( 2*t - 3) + 1; // H00 = 2t^3 - 3t^2+1 Hs[j,1] := t_t*( t - 2) + t; // H10 = t^3 - 2t^2+t Hs[j,2] := t_t*(-2*t + 3); // H01 =-2t^3 + 3t^2 Hs[j,3] := t_t*( t - 1); // H11 = t^3 - t^2 end; // Формируем массив результата if Closed then begin L := N; SetLength(Result, L*PointsPerSegment); end else begin L := N - 1; SetLength(Result, L*PointsPerSegment+1); end; Idx := 0; // Генерация сплайна с использованием кэша for i := 0 to L-1 do begin P1 := Points[i]; // Соседние точки для расчета касательных if Closed then begin P0 := Points[(i-1+N) mod N]; P2 := Points[(i+1) mod N]; P3 := Points[(i+2) mod N]; end else begin if i=0 then begin P0 := Points[0]; P2 := Points[1]; end else P0 := Points[i-1]; if i=L-1 then begin P2 := Points[N-1]; P3 := Points[N-1]; end else begin P2 := Points[i+1]; P3 := Points[i+2]; end; end; // Касательные векторы с учетом натяжения M0 := (P2 - P0)*Tension*0.5; M1 := (P3 - P1)*Tension*0.5; // Генерация точек сегмента for j := 0 to PointsPerSegment-1 do begin if not Closed and (i=L-1) and (j=PointsPerSegment) then break; // Для последней точки разомкнутой кривой // Вычисление точки кривой по функции Cubic Hermite Splines pt.X := Hs[j,0]*P1.X + Hs[j,1]*M0.X + Hs[j,2]*P2.X + Hs[j,3]*M1.X; pt.Y := Hs[j,0]*P1.Y + Hs[j,1]*M0.Y + Hs[j,2]*P2.Y + Hs[j,3]*M1.Y; Result[Idx] := Pt; Inc(Idx); end; end; // Добавление последней точки для разомкнутой кривой if not Closed then Result[Idx] := Points[N-1]; end; |
Сплайн Эрмита: Онлайн
Предпринял попытку показать, как строится сплайн Эрмита в динамике. Наверное, получилось всё равно не слишком информативно, поэтому буду рад всем пожеланиям для улучшения наглядности. Я за популяризацию математики!
Тут рассмотрен «хороший» случай, когда у нас есть интервал (P1-P2), у которого есть и предыдущая точка (P0) и следующая (P2). То есть, дело происходит где-то посреди гипотетической ломаной. Необходимо выставить параметры натяжения и числа точек внутри отрезка и нажать кнопку Построить сплайн.
Картинка становится интересней, когда натяжение велико (>6) и число точек, скажем, 200.
За красные точки можно таскать!
Демонстрационная прога
Добавить точку можно двойным кликом. Убрать точку можно двойным кликом по существующей точке. Предупреждений про удаление нет, откатить обратно нельзя, так что, аккуратней. Это ж демка, блин, а не графический редактор!
При наведении мышкой на точку измениться курсор. Это значит, точку можно таскать. Не надо утаскивать за пределы экрана, масштаба тут нет, поэтому добраться до неё может оказаться невозможным.
Если были произведены изменения, при закрытии прога спросит, сохранить ли точки. Если они будут сохранены, при следующем открытии сохранённые точки появятся на экране. Они хранятся в файлике points.sav рядом с прогой, можно его удалить и точки при открытии будут в состоянии по умолчанию.
| Tension | Параметр натяжения кривой. В реальности передаётся значение, делённое на 10. В идеале должен быть в интервале 0..1 (0..10 в проге), 0-форма не меняется, 1-максимальное натяжение. Но можно задавать и больше, ради эксперимента. |
| Points | Количество точек на сегмент. |
| Squared | Если флаг установлен, координаты рассчитываются исходя из того, что масштаб по X и Y одинаков. Если снять флаг, изображение вытянется по горизонтали (или вертикали). |
| Closed | Если флаг установлен, значит имеем дело с замкнутой кривой. Полезно отследить поведение кривой в замкнутом и разомкнутом состоянии |
| Draw Curve | Нарисовать кривую силами GDI+ на тех же исходных точках (DrawCurve). Для сравнения. DrawCurve, оказывается, не использует кубические сплайны Эрмита. |
| Line Only | Рисовать сглаженную кривую только линиями, без точек. Полезно, чтобы оценить форму кривой. |
Скачать
Друзья, спасибо за внимание!
Исходник (zip) 69 Кб. Delphi XE 7
Исполняемый файл (zip) 933 Kб (Скомпилирован в XE 7)
PS. Ссылки по теме
Математический сплайн (wiki)
Сплайн_Эрмита (wiki)
Моделирование кривых и поверхностей. Сплайны (pdf)
Кубический сплайн (pdf)
Преобразование кубического эрмитова представления кривых и поверхностей в кубический B-сплайн (pdf)
Роджерс Д., Адамс Дж. Математические основы машинной графики (pdf) 5.74 Mb
On Cardinal Spline Interpolation (pdf) 623 Kb
КУРС «ЧИСЛЕННЫЕ МЕТОДЫ» Интерполяция кубическими сплайнами (pdf) 559 Kb
Построение сплайнов с использованием OpenGL (pdf) 2.39 Mb