Всем привет!

Как и понятное по заголовку, это продолжение предыдущей статьи с той же темой, но теперь мы будем рассматривать аппроксимационные кривые.

Основная часть

В прошлой статье я рассматривал семейство интерполяционных кривых вида:

$y(x)=\frac{\sum_{i=1}^{n}y_{i}\cdot k_{i}}{\sum_{i=1}^{n}k_{i}}.$

Функции параметра для интерполяции принимают значения от до и возвращают значения от до . Как я и обещал в предыдущей статье, противоположные функции будут в этой статье и текущая тема будет уже не интерполяция, а аппроксимация.

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

А теперь, барабанная дробь, это семейство не только аппроксимирует, но и позволяет регулировать точность аппроксимации! Это легко доказывается. При $\textstyle x\in\left[0;1\right], k\geq0,$ справедливо неравенство $x^{n+k}\leq x^{n}$ , то есть график первой функции сильнее прижат к оси . Если провести аналогию с расстоянием ( принимаем за расстояние до точки), то на первом графике расстояние влияет сильнее, но одинаковое в и . Но важно помнить, что при очень высокой степени точность на самих точках будет высокая, но поблизости значение функции будет прижато к значению точки.

Настало время перейти непосредственно к рассмотрению функций параметра. Функции могут быть почти те же, что и для интерполяции, только будет отличаться их применение — нужно будет вычитать их значение из максимального. Примерами таких функций являются, аналогично с интерполяцией, $x^{n}$ , $\sin\left(x\right)$ , $n^{x}$ , $\tan\left(x\right)$ и другие.

Квадрат

Не отступая от плана предыдущей статьи, первая функция на очереди — степенная. Если для интерполяции мы делили, то здесь нам придется умножать. Сам вид тоже чуточку меняется и, соответственно, выглядит так.

$k=\left(1-\left(\frac{x-x_{i}}{x_{n}-x_{1}}\right)^2\right)^{p}\implies k=\left(\left(x_{n}-x_{1}\right)^{2}-\left(x-x_{i}\right)^{2}\right)^{p},$

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

public static float power(Vector2[] points, float x, float n, float p)
{
  float b = 0, y = 0;
  for (int i = 0; i < points.Length; i++)
  {
    float k = Mathf.Pow(Mathf.Pow(Mathf.Abs(points[points.Length - 1].x - points[0].x), n) - Mathf.Pow(Mathf.Abs(x - points[i].x), n), p);
    float y += a[i].y * k;
    float b += k;
  }
  return y / b;
}

Синус/Косинус

косинус, растяжение 1, сдвиг 1, точность 20

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

Итак, наша функция должна возвращать при входном , и наоборот, соответственно. Параметров всего два — растяжение и сдвиг. Максимальный полупериод может быть равен , поэтому формула растяжения будет такой: $\cos\left(k\cdot\pi\cdot x\right)$ , где — коэффициент $\in\left(0;1\right]$ . С добавлением сдвига она превращается в $\cos\left(k_{1}\cdot\pi\cdot x-k_{2}\cdot\pi\right)$ . Теперь нам надо привести формулу к нужному виду, поэтому окончательный вид немного другой:

$k=\left(\frac{\cos\left(k_{1}\cdot\pi\cdot x-k_{2}\cdot\pi\right)-\cos\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)}{\cos\left(k_{2}\cdot\pi\right)-\cos\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)}\right)^{p}\implies$ $\implies k=\left(\cos\left(k_{1}\cdot\pi\cdot x-k_{2}\cdot\pi\right)-\cos\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)\right)^{p}.$

Коэффициенты ограничены — $k_{1}\in\left(0;1\right]$ , а $k_{2}\in\left[0;1\right]/\left\{\frac{k_{1}}{2}\right\}$ . Остается написать код для синусоидой аппроксимации.

public static float sinus(Vector2[] points, float x, float k1, float k2, float p)
{
  float b = 0, y = 0;
  for (int i = 0; i < points.Length; i++)
  {
    float k = Mathf.Pow(Mathf.Cos(k1 * Mathf.PI * Mathf.Abs(x - points[i].x) / (points[points.Length - 1].x - points[0].x) - k2 * Mathf.PI) - Mathf.Cos(k1 * Mathf.PI - k2 * Mathf.PI), p);
    y += points[i].y * k;
    b += k;
  }
  return y / b;
}

Экспонента

экспонента, второй коэффициент 4, точность 30

Настало время рассмотреть функции вида $n^{x}$ . В рамках данной статьи, хочу заметить, что они не очень удобные, поскольку при некоторых значениях параметра контролировать степень аппроксимации почти невозможно. Но для программирования это может оказаться полезным, поскольку возводить в очень большие степени не потребуется. Для функции такого вида есть несколько вариантов, подходящих для аппроксимации, но я выберу лучшие два. Сверху у нас должна быть степень наибольшая при $\Delta x\to0$ и наименьшая при $\Delta x\to\max$ . Для этого можно использовать квадратичную функцию, с которой мы уже познакомились в этой статье. По свойству степеней у нас степени перемножаются:

$k=c^{\displaystyle k_{квадрат}\cdot p}.$

Если раскрыть, получается три коэффициента:

$k=k_{1}^{\displaystyle k_{2}\cdot\left(1-\left(\frac{x-x_{i}}{x_{n}-x_{1}}\right)^{2}\right)^{p}}.$

Второй коэффициент не стоит делать слишком большим, иначе будет трудно контролировать аппроксимацию. В плане кода это будет два возведения в степень, но не очень большие. Ну и последний шаг для этого блока — код.

public static float exponent(Vector2[] points, float k1, float k2, float p)
{
  float b = 0, y = 0;
  for (int i = 0; i < points.Length; i++)
  {
    float k = Mathf.Pow(k1, k2 * Mathf.Pow(1 - (x - points[i].x) * (x - points[i].x) / (points[points.Length - 1].x - points[0].x) / (points[points.Length - 1].x - points[0].x), p));
    y += points[i].y * k;
    b += k;
  }
  return y / b;
}

Тангенс/Котангенс

тангенс, растяжение 0.1, сдвиг 0, точность 18

Продолжаем тригонометрию. Из основных тригонометрических функций осталась одна — тангенс. Для тангенса все аналогично с синусом, вид такой:

$k=\left(\frac{\tan\left(k_{1}\cdot\pi\cdot x-k_{2}\cdot\pi\right)-\tan\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)}{-\tan\left(k_{2}\cdot\pi\right)-\tan\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)}\right)^{p}\implies$ $\implies k=\left(\tan\left(k_{1}\cdot\pi\cdot x-k_{2}\cdot\pi\right)-\tan\left(k_{1}\cdot\pi-k_{2}\cdot\pi\right)\right)^{p}$

Ограничения коэффициентов — $k_{1}\in\left(0;\frac{1}{2}\right]$ , $k_{2}\in\left[0;\frac{1}{2}\right)$ , $k_{1}-k_{2}\ne \frac{n\cdot\pi}{2}$ .

Для нужной функции в коде меняем только косинус на тангенс в коде из блока про синус. Да уж, халявный блок получается, по готовенькому прошли).

public static float tangent(Vector2[] points, float x, float k1, float k2, float p)
{
  float b = 0, y = 0;
  for (int i = 0; i < points.Length; i++)
  {
    float k = Mathf.Pow(Mathf.Tan(k1 * Mathf.PI * Mathf.Abs(x - points[i].x) / (points[points.Length - 1].x - points[0].x) - k2 * Mathf.PI) - Mathf.Tan(k1 * Mathf.PI - k2 * Mathf.PI), p);
    y += points[i].y * k;
    b += k;
  }
  return y / b;
}

Заключение

Ну и куда же без кастомной функции и самой маленькой оптимизации кода:

public struct Approximation
{
  public static float power(Vector2[] points, float x, float n, float p)
  {
    float b = 0, y = 0;
    float m = Mathf.Pow(Mathf.Abs(points[points.Length - 1].x - points[0].x), n);
    for (int i = 0; i < points.Length; i++)
    {
      float k = Mathf.Pow(m - Mathf.Pow(Mathf.Abs(x - points[i].x), n), p);
      float y += a[i].y * k;
      float b += k;
    }
    return y / b;
  }
  
  public static float sinus(Vector2[] points, float x, float k1, float k2, float p)
  {
    float b = 0, y = 0;
    float c = Mathf.Cos(k1 * Mathf.PI - k2 * Mathf.PI), m = points[points.Length - 1].x - points[0].x;
    for (int i = 0; i < points.Length; i++)
    {
      float k = Mathf.Pow(Mathf.Cos(k1 * Mathf.PI * Mathf.Abs(x - points[i].x) / m - k2 * Mathf.PI) - c, p);
      y += points[i].y * k;
      b += k;
    }
    return y / b;
  }
  
  public static float exponent(Vector2[] points, float k1, float k2, float p)
  {
    float b = 0, y = 0;
    float m = (points[points.Length - 1].x - points[0].x) * (points[points.Length - 1].x - points[0].x);
    for (int i = 0; i < points.Length; i++)
    {
      float k = Mathf.Pow(k1, k2 * Mathf.Pow(1 - (x - points[i].x) * (x - points[i].x) / m, p));
      y += points[i].y * k;
      b += k;
    }
    return y / b;
  }
  
  public static float tangent(Vector2[] points, float x, float k1, float k2, float p)
  {
    float b = 0, y = 0;
    float t = Mathf.Tan(k1 * Mathf.PI - k2 * Mathf.PI), m = points[points.Length - 1].x - points[0].x;
    for (int i = 0; i < points.Length; i++)
    {
      float k = Mathf.Pow(Mathf.Tan(k1 * Mathf.PI * Mathf.Abs(x - points[i].x) / m - k2 * Mathf.PI) - t, p);
      y += points[i].y * k;
      b += k;
    }
    return y / b;
  }
  
  public delegate float function(float f);
  
  public static delegat(Vector2[] points, float x, function func)
  {
    float b = 0, y = 0;
    float m = points[points.Length - 1].x - points[1].x;
    for (int i = 0; i < points.Length; i++)
    {
      float k = f(Mathf.Abs(points[i].x - x) / m);
      y += points[i].y * k;
      b += k;
    }
    return y / b;
  }
}

В данной статье рассмотрены лишь немногие функции, остальные остаются можете рассмотреть самостоятельно. Код в главной части статьи приведен полностью по материалу, то есть нет оптимизаций, из-за которых код будет по виду отличаться от самой формулы. Используется класс для точек от движка Unity, вещественные числа типа float именно из-за движка и применяются, для чистого C# нужно использовать double.

Хочу еще немножко сказать про аппроксимацию. При работе с большим количеством данных, важно помнить, что количество итераций растет линейно. Данное семейство также позволяет экспериментировать с различными функциями — забавно посмотреть на аппроксимацию с функциями вида $f(x)=1-2\cdot x^{4}+x^{2}$ .

Возможностей для применения такой аппроксимации много, но три я хочу выделить. Первая, и самая очевидная, — это аппроксимация движения между точками. Функционал позволяет настроить это так, как вам нужно. Вторая — это аппроксимация данных, предсказание новых по старым. На мой взгляд, это самая перспективная область для аппроксимации. Ну и третья — обработка изображений, функция сглаживания. При мощном компе можно обрабатывать большие изображения, либо ограничиваться ближайшими точками. Можно увеличивать изображение, но тогда нельзя забывать оставлять узловые значения, чтобы оно не теряло четкость без надобности.

Надеюсь, вам где-нибудь это пригодится)

Кривые и что это такое ч.3