Кубические кривые Безье - получить Y для заданного X - особый случай, когда X контрольных точек увеличивается

Question

Я прочитал несколько обсуждений относительно нахождения Y в X для кубической кривой Безье, а также прочитал эту статью по этому вопросу.

Мой случай более ограничен, чем общий, и мне интересно, есть ли лучшее решение, чем общие, упомянутые в приведенных выше обсуждениях.

Мое дело:

• X Значение различных контрольных точек увеличивается. То есть: X3 > X2 > X1 > X0,
• Кроме того, в результате вышеизложенного, X(t) также строго монотонно возрастает.

Есть ли эффективный алгоритм, который учитывает такие ограничения?

Answer 1

Прежде всего: этот ответ работает только потому, что ваше ограничение контрольной точки означает, что мы всегда имеем дело с параметрическим эквивалентом нормальной функции. Это, очевидно, то, что вам нужно в этом случае, но любой, кто найдет этот ответ в будущем, должен знать, что этот ответ основан на предположении, что для любого данного значения x существует только одно значение y. Что, конечно, абсолютно не верно для кривых Безье в целом.

С учетом сказанного мы знаем, что, хотя мы выразили эту кривую как параметрическую кривую в двух измерениях, мы имеем дело с кривой, которая для всех намерений и целей должна иметь некоторую неизвестную функцию вида y = f(x), Мы также знаем, что до тех пор, пока мы знаем значение "t", которое однозначно принадлежит конкретному x (что имеет место только из-за вашего строго монотонно возрастающего свойства коэффициентов), мы можем вычислить y как y = By(t) поэтому возникает вопрос: можем ли мы вычислить t значение, которое нам нужно подключить к By(t), учитывая некоторые известные x значение?

На что ответ: да, мы можем.

Во-первых, любой x Можно сказать, что значение, с которого мы начинаем x = Bx(t) так, учитывая, что мы знаем x, мы должны быть в состоянии найти соответствующее значение (я) t что приводит к этому x,

давайте посмотрим на функцию для x(t):

x(t) = a(1-t)³ + 3b(1-t)²t + 3c(1-t)t² + dt³

Мы можем переписать это в простую полиномиальную форму как:

x(t) = (-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + a

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

(-a + 3b- 3c + d)t³ + (3a - 6b + 3c)t² + (-3a + 3b)t + (a-x) = 0

Вы можете спросить: "Куда делись все -x для других значений a, b, c и d?" и ответ здесь заключается в том, что все они отменяются, поэтому единственное, что нам на самом деле нужно вычесть, это тот, который полностью в конце. Handy!

Итак, теперь просто... решить это уравнение: мы знаем все, кроме t Нам просто нужно некоторое математическое понимание, чтобы сказать нам, как это сделать.

... Конечно, "просто" не является правильным определителем здесь, нет ничего "просто" в поиске корней кубической функции, но, к счастью, Джеролано Кардано заложил основы для определения корней еще в 16 веке, используя комплексные числа. Еще до того, как кто-то изобрел комплексные числа. Совершенно подвиг! Но это программный ответ, а не урок истории, так что давайте приступим к реализации:

Учитывая некоторое известное значение для x и знание наших координат a, b, c и d, мы можем реализовать наш поиск корня следующим образом:

// Find the roots for a cubic polynomial with bernstein coefficients
// {pa, pb, pc, pd}. The function will first convert those to the
// standard polynomial coefficients, and then run through Cardano's
// formula for finding the roots of a depressed cubic curve.
double[] findRoots(double x, double pa, double pb, double pc, double pd) {
  double
    pa3 = 3 * pa,
    pb3 = 3 * pb,
    pc3 = 3 * pc,
    a = -pa  +   pb3 - pc3 + pd,     
    b =  pa3 - 2*pb3 + pc3, 
    c = -pa3 +   pb3, 
    d =  pa  -     x;

  // Fun fact: any Bezier curve may (accidentally or on purpose)
  // perfectly model any lower order curve, so we want to test 
  // for that: lower order curves are much easier to root-find.
  if (approximately(a, 0)) {
    // this is not a cubic curve.
    if (approximately(b, 0)) {
      // in fact, this is not a quadratic curve either.
      if (approximately(c, 0)) {
        // in fact in fact, there are no solutions.
        return new double[]{};
      }
      // linear solution:
      return new double[]{-d / c};
    }
    // quadratic solution:
    double
      q = sqrt(c * c - 4 * b * d), 
      b2 = 2 * b;
    return new double[]{
      (q - c) / b2, 
      (-c - q) / b2
    };
  }

  // At this point, we know we need a cubic solution,
  // and the above a/b/c/d values were technically
  // a pre-optimized set because a might be zero and
  // that would cause the following divisions to error.

  b /= a;
  c /= a;
  d /= a;

  double
    b3 = b / 3,
    p = (3 * c - b*b) / 3, 
    p3 = p / 3, 
    q = (2 * b*b*b - 9 * b * c + 27 * d) / 27, 
    q2 = q / 2, 
    discriminant = q2*q2 + p3*p3*p3, 
    u1, v1;

  // case 1: three real roots, but finding them involves complex
  // maths. Since we don't have a complex data type, we use trig
  // instead, because complex numbers have nice geometric properties.
  if (discriminant < 0) {
    double
      mp3 = -p/3,
      r = sqrt(mp3*mp3*mp3), 
      t = -q / (2 * r), 
      cosphi = t < -1 ? -1 : t > 1 ? 1 : t, 
      phi = acos(cosphi), 
      crtr = crt(r), 
      t1 = 2 * crtr;
    return new double[]{
      t1 * cos(phi / 3) - b3,
      t1 * cos((phi + TAU) / 3) - b3,
      t1 * cos((phi + 2 * TAU) / 3) - b3
    };
  }

  // case 2: three real roots, but two form a "double root",
  // and so will have the same resultant value. We only need
  // to return two values in this case.
  else if (discriminant == 0) {
    u1 = q2 < 0 ? crt(-q2) : -crt(q2);
    return new double[]{
      2 * u1 - b3,
      -u1 - b3
    };
  }

  // case 3: one real root, 2 complex roots. We don't care about
  // complex results so we just ignore those and directly compute
  // that single real root.
  else {
    double sd = sqrt(discriminant);
    u1 = crt(-q2 + sd);
    v1 = crt(q2 + sd);
    return new double[]{u1 - v1 - b3};
  }
}

Хорошо, это довольно кусок кода, с довольно многими дополнениями:

• crt() это функция cuberoot Мы на самом деле не заботимся о комплексных числах в этом случае, так что проще всего реализовать это с помощью def, или макроса, или троичного, или любого другого сокращения, которое предлагает ваш язык по вашему выбору: crt(x) = x < 0 ? -pow(-x, 1f/3f) : pow(x, 1f/3f);,
• tau просто 2π. Это полезно, когда вы занимаетесь программированием геометрии.
• approximately это функция, которая сравнивает значение с очень маленьким интервалом вокруг цели, потому что числа с плавающей запятой IEEE являются рывками. В основном мы говорим о approximately(a,b) = return abs(a-b) < 0.000001 или что-то.

Остальное должно быть довольно самоочевидным, если немного Java-esque (я использую Processing для такого рода вещей).

С помощью этой реализации мы можем написать нашу реализацию, чтобы найти y по заданному x. Что немного сложнее, чем просто вызов функции, потому что кубические корни - сложные вещи. Вы можете получить до трех корней назад. Но только один из них будет применим к нашему "t-интервалу" [0,1]:

double x = some value we know!
double[] roots = getTforX(x);
double t;
if (roots.length > 0) {
  for (double _t: roots) {
    if (_t<0 || _t>1) continue;
    t = _t;
    break;
  }
}

И все, мы закончили: теперь у нас есть значение "t", которое мы можем использовать, чтобы получить соответствующее значение "y".

Answer 2

Если бинарный поиск слишком сложен, все еще есть O(1) подход, но его довольно ограниченный. Я полагаю, вы используете 4 контрольной точки (p0(x0,y0),p1(x1,y1),p2(x2,y2),p3(x3,y3)б) кубический безье параметризован t в промежутке [0.0 , 1.0] так:

t = 0.0 -> x(t) = x0, y(t) = y0;
t = 1.0 -> x(t) = x3, y(t) = y3;

Сначала давайте на некоторое время забудем о Безье и вместо этого используем кривую Кэт-Ромла, которая является просто альтернативным способом представления этой же кривой. Для преобразования между 2 кубиками используйте эти:

// BEzier to Catmull-Rom
const double m=6.0;
X0 = x3+(x0-x1)*m; Y0 = y3+(y0-y1)*m;
X1 = x0;           Y1 = y0;
X2 = x3;           Y2 = y3;
X3 = x0+(x3-x2)*m; Y3 = y0+(y3-y2)*m;

// Catmull-Rom to Bezier
const double m=1.0/6.0;
x0 = X1;           y0 = Y1;
x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
x3 = X2;           y3 = Y2;

где (xi,yi) контрольные точки Безье и (Xi,Yi) являются точками Кэтмулла-Рома. Теперь, если X Расстояние между всеми контрольными точками имеет одинаковое расстояние:

(X3-X2) == (X2-X1) == (X1-X0)

тогда X координата является линейной с t, Это означает, что мы можем вычислить t прямо из X:

t = (X-X1)/(X2-X1);

Теперь мы можем вычислить Y для любого X непосредственно. Поэтому, если вы можете выбрать контрольные точки, выберите их так, чтобы они соответствовали условию X-расстояния.

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