Доброго времени суток, харбачитатель.
Так начинается статья, которая представляет сообществу первый, опубликованный здесь, алгоритм интерполяции:
В этой статье мне хотелось бы рассказать об одном придуманном алгоритме (или скорее всего — переизобретённом велосипеде) построения плавного графика по заданным точкам, используя кривые Безье. Статья была написана под влиянием вот этой статьи...
Дело в том, что автор этой статьи пошел не прямым путем, а окольным, используя громоздкую тригонометрию половинных углов для вычисления тангенса угла φ. Кроме того автор не уделил внимание анализу длин опорных отрезков, просто приняв, что данные равномерно распределены по оси Х с фиксированным шагом, и поэтому длины опорных отрезков тоже можно взять фиксированными по оси X.
На самом деле, для того, чтобы вычислить координаты опорных точек не нужна тригонометрия, не нужно вычислять угол φ и даже не нужен коэффициент k. Всё что нужно — это знание теоремы Пифагора и немного смекалки в определении пропорций между отрезками. При этом у меня получился алгоритм, не требующий равномерности шага и даже строгой направленности вдоль оси X. Вы можете направлять данные произвольно, хоть вдоль оси Y, хоть по кругу.
Для понимания алгоритма вычисления координат опорных точек построим рисунок, аналогичный соответствующему рисунку из статьи первой версии.
Я только перенес точку D на отрезок AB таким образом, чтобы отрезок CD был также наклонён к оси X, как и отрезок B1C1. Нетрудно понять, что при этом треугольник DAC является равнобедренным. Также заметим, что длина отрезка CD равна двойной величине по длине опорного отрезка справа от точки A, при котором интерполяционная дуга будет близка к окружности. На рисунке выше такая длина обозначена двумя серыми линиями, перпендикулярными красной линии с опорными отрезками и отрезку CD. Примем длину CD в 2 относительные единицы. Также сразу заметим, что максимально возможная длина опорного отрезка слева в таких же единицах будет равна отношению |AB|/|AC|. Соответственно в случае, когда |AB|=|AC|, максимальная длина опорного отрезка слева тоже будет равна 1. Обозначим через m требуемую для построения относительную длину опорных отрезков. Обычно при использовании кубических кривых Безье m=0.333, а для квадратичных m=0.5.
Теперь мы можем написать формулы для расчета координат опорных точек слева и справа:
XB1 = XA-m*ΔXCD*|AB|/|AC|/2 | (1) | |
YB1 = YA-m*ΔYCD*|AB|/|AC|/2 | (2) | |
XC1 = XA+m*ΔXCD/2 | (3) | |
YC1 = YA+m*ΔYCD/2 | (4) |
Здесь
X и Y — это координаты, соответствующие обозначениям, точек;
ΔXCD=XС-XD и ΔYCD=YС-YD (нет модуля, значения могут быть больше или меньше нуля).
Абсолютные значения модулей длин отрезков легко определяются с помощью теоремы Пифагора. Поэтому по сути нам здесь может представлять некоторую сложность определение координат точки D.
Давайте примем всю длину отрезка AB, на которой располагается точка D, за 1. Тогда длина отрезка AD относительно AB (обозначим через d) d=|AC|/|AB| и координаты точки D:
XD = XA-d*ΔXAB | (5) | |
YD = YA-d*ΔYAB | (6) |
И несколько оптимизируем формулы (1), (2), (3) и (4) сделав замену переменных, используя d и m2=m/2:
XB1 = XA-m2*ΔXCD /d | (7) | |
YB1 = YA-m2*ΔYCD /d | (8) | |
XC1 = XA+m2*ΔXCD | (9) | |
YC1 = YA+m2*ΔYCD | (10) |
И последний штрих, который на мой взгляд надо добавить в алгоритм — это устранение завалов, которые появляются при сильных углах наклона опорных отрезков, в том числе при большой разнице масштаба по оси Y по сравнению с масштабом по оси X. В такой ситуации кривая может через чур сильно начать вытягиваться направо или налево. Разумно полагать, что всё-таки максимумы и минимумы в природе более-менее симметричны и справа и слева. Поэтому мы можем ввести более мягкое условие, чем требовать от исходных данных равномерный шаг. К тому же интерполяция зачастую и необходима для того, чтобы привести исходные данные с неравномерным шагом к равномерному шагу. Будет логично длину опорных отрезков справа и слева сделать одинаковыми, ну или близкими по величине. Для начала мы сделаем замену переменных в формулах 7 и 8, обозначив через m1=m2/d:
XB1 = XA-m1*ΔXCD | (11) | |
YB1 = YA-m1*ΔYCD | (12) |
Рассчитав m1 и m2, мы проверяем, чтобы большее из них не было слишком велико. Можно установить требование, скажем чтобы больший из опорных отрезков был не более, чем на 7% больше меньшего. Визуально разница длин опорных отрезков на 5% совсем не заметна, при разнице в 10% уже немного заметна асимметрия.
Пример текста функции расчета опорных точек на MS VFP
function BezierOpt(Xb,Yb,Xa,Ya,Xc,Yc,m1,m2)
** Определение опорных точек. Версия 2
** A - средняя точка. Теория здесь https://habr.com/ru/sandbox/221744/
** Determining the reference points
** A - the middle point.
** Determinacion de los puntos de referencia
** A - el punto medio.
Xab=m.Xa-m.Xb
Yab=m.Ya-m.Yb
Xca=m.Xc-m.Xa
Yca=m.Yc-m.Ya
AB=sqrt(m.Xab^2+m.Yab^2)
AC=sqrt(m.Xca^2+m.Yca^2)
t=iif(empt(m.AB),m.AB,m.AC/m.AB)
Xd=m.Xa-m.Xab*m.t
Yd=m.Ya-m.Yab*m.t
Xcd=m.Xc-m.Xd
Ycd=m.Yc-m.Yd
t1=iif(empt(m.t),m.t,m.m1/m.t/2)
t2=m.m2/2
Xb1=m.Xa-m.t1*m.Xcd
Yb1=m.Ya-m.t1*m.Ycd
Xc1=m.Xa+m.t2*m.Xcd
Yc1=m.Ya+m.t2*m.YcdСобственно это и было моим ответом автору первой версии этого алгоритма. Работу алгоритма в режиме онлайн вы можете проверить на сайте.
И стоит отметить, что автор первой версии абсолютно правильно заметил, что для того, чтобы по оси X интерполируемые значения не направлялись в обратную сторону, следует ограничивать по длине опорные отрезки, так, чтобы проекция опорного отрезка на ось X на каждом интервале не превышала 0.5 ΔX.
А как с построением опорных отрезков у офисного пакета американской корпорации Microsoft или версия алгоритма 2.1
Слабое место алгоритма версии 2, в том, что приходится ограничивать длины слишком длинных опорных отрезков. Зачем изобретать велосипед? — подумал я и решил посмотреть, а как строит опорные отрезки при рисовании кривых известный офисный пакет. Оказалось, что наклон опорных отрезков у Microsoft совсем не такой, как предлагает автор первой версии. Наклон опорных отрезков Microsoft совпадает с наклоном основания треугольника ABC. И это правильно — не надо делить вершину A пополам, а потом думать, как ограничить получившийся неадекватно длинный опорный отрезок. При таком наклоне, как у Microsoft всё получается гармонично. Но Microsoft считает, что длины отрезков слева и справа должны быть одинаковыми... Интуитивно длины опорных отрезков всё-таки должны зависеть от левой и правой сторон треугольника. Если провести высоту к основанию треугольника, то такая зависимость будет весьма сильной, возможно излишне сильной. Мы уже сталкивались с завалами, а потом добавляли в алгоритм ограничение встретившегося слишком длинного опорного отрезка. Поэтому на этот раз я пойду на компромисс и проведу биссектрису из вершины A. При этом зависимость от сторон треугольника остается, но она более слабая, чем если бы мы проводили перпендикулярную к основанию высоту.
Координаты точки D:
XD=XB+(XC-XB)*t | (13) | |
YD=YB+(YC-YB)*t | (14) |
где t — это относительная длина отрезка |AB| к сумме |AB|+|AC|, t=|AB|/(|AB|+|AC|).
И координаты концов опорных отрезков:
XB1 = XA-m*ΔXDB | (15) | |
YB1 = YA-m*ΔYDB | (16) | |
XC1 = XA+m*ΔXCD | (17) | |
YC1 = YA+m*ΔYCD | (18) |
Пример текста функции расчета опорных точек на MS VFP
function BezierOpt(Xb,Yb,Xa,Ya,Xc,Yc,m1,m2)
** Определение опорных точек. Версия 2.1
** A - средняя точка. Теория здесь https://habr.com/ru/articles/831662/
** Determining the reference points
** A - the middle point.
** Determinacion de los puntos de referencia
** A - el punto medio.
Xab=m.Xa-m.Xb
Yab=m.Ya-m.Yb
Xca=m.Xc-m.Xa
Yca=m.Yc-m.Ya
AB=sqrt(m.Xab^2+m.Yab^2)
AC=sqrt(m.Xca^2+m.Yca^2)
t=iif(empt(m.AB+m.AC),m.AB,m.AB/(m.AB+m.AC))
Xd=m.Xb+(m.Xc-m.Xb)*m.t
Yd=m.Yb+(m.Yc-m.Yb)*m.t
Xb1=m.Xa+m.m1*(m.Xb-m.Xd)
Yb1=m.Ya+m.m1*(m.Yb-m.Yd)
Xc1=m.Xa+m.m2*(m.Xc-m.Xd)
Yc1=m.Ya+m.m2*(m.Yc-m.Yd)Работу алгоритма в режиме онлайн вы можете проверить у меня на сайте. Кстати сайт работает на https.net сервере собственной разработки под Windows.
От двумерной к одномерной интерполяции.
Если стоит задача интерполировать данные однонаправленной функции или процесса, т.е. функции зависящей от времени, то чтобы интерполируемые значения гарантированно не возвращались назад по оси X для быстро растущих или падающих зависимостей, необходимо установить ограничение на длины опорных отрезков, таким образом, чтобы проекция этих отрезков на ось X не превышала половины шага исходных данных по оси X, если такое будет иметь место.
Я попробовал взять за основу для модели одномерной интерполяции обе версии 2 и 2.1 и пришел к выводу, что версия 2.1 всё-таки дает более плавные однонаправленные кривые, чем 2. Поэтому далее рассмотрим получение формулы требуемого ограничения для версии 2.1.
Собственно само ограничение вытекает из формул (15) и (17):
m*|ΔXDB| ≤ 0.5*|ΔXAB| | (19) | |
m*|ΔXCD| ≤ 0.5*|ΔXAC| | (20) |
или
mB = |ΔXAB/ΔXDB|/2 | (21) | |
mC = |ΔXAC/ΔXCD|/2 | (22) |
где mB и mC максимальные значения, ограничивающие коэффициент m слева и справа.
Пример текста функции расчета опорных точек на MS VFP
function BezierOpt(Xb,Yb,Xa,Ya,Xc,Yc,m1,m2)
** Определение опорных точек для одномерной интерполяции
** A - средняя точка. Теория здесь https://habr.com/ru/articles/831662/
** Determining the reference points
** A - the middle point.
** Determinacion de los puntos de referencia
** A - el punto medio.
Xab=m.Xa-m.Xb
Yab=m.Ya-m.Yb
Xca=m.Xc-m.Xa
Yca=m.Yc-m.Ya
AB=sqrt(m.Xab^2+m.Yab^2)
AC=sqrt(m.Xca^2+m.Yca^2)
t=iif(empt(m.AB+m.AC),m.AB,m.AB/(m.AB+m.AC))
Xd=m.Xb+(m.Xc-m.Xb)*m.t
Yd=m.Yb+(m.Yc-m.Yb)*m.t
Xbd=m.Xb-m.Xd
Xcd=m.Xc-m.Xd
mb=abs(m.Xab/m.Xbd)/2
if m.mb>m.m1
mb=m.m1
endi
mc=abs(m.Xca/m.Xcd)/2
if m.mc>m.m2
mc=m.m2
endi
Xb1=m.Xa+m.mb*m.Xbd
Yb1=m.Ya+m.mb*(m.Yb-m.Yd)
Xc1=m.Xa+m.mc*m.Xcd
Yc1=m.Ya+m.mc*(m.Yc-m.Yd)Работу алгоритма для интерполяции одномерных однонаправленных функций можно проверить здесь.
