Как вычислить вершину параболы с учетом трех точек
У меня есть три точки X/Y, которые образуют параболу. Мне просто нужно вычислить, какая вершина параболы проходит через эти три точки. Желательно, чтобы это был быстрый способ, поскольку мне нужно сделать МНОГО таких вычислений!
Веб-сайт "Спроси ученого" дает такой ответ:
Общий вид параболы задается уравнением: A * x^2 + B * x + C = y, где A, B и C - произвольные действительные постоянные. У вас есть три пары точек, которые являются (x,y) упорядоченными парами. Подставим значения x и y каждой точки в уравнение для параболы. Вы получите три ЛИНЕЙНЫХ уравнения с тремя неизвестными, тремя константами. Затем вы можете легко решить эту систему из трех уравнений для значений A, B и C, и вы получите уравнение параболы, которая пересекает ваши 3 точки. В вершине, где первая производная равна 0, маленькая алгебра дает: ( -B/2A, C - B^2/4A) для вершины.
Было бы неплохо увидеть реальный код, который выполняет эти вычисления в C# или C++. Кто-нибудь?
9 ответов
На самом деле это простая задача линейной алгебры, поэтому вы можете выполнять вычисления символически. Подставив значения x и y трех ваших точек, вы получите три линейных уравнения с тремя неизвестными.
A x1^2 + B x1 + C = y1
A x2^2 + B x2 + C = y2
A x3^2 + B x3 + C = y3
Простой способ решить эту проблему - инвертировать матрицу.
x1^2 x1 1
x2^2 x2 1
x3^2 x3 1
и умножить его на вектор
y1
y2
y3
В результате... ладно, не все так просто;-) Я сделал это в Mathematica, а вот формулы в псевдокоде:
denom = (x1 - x2)(x1 - x3)(x2 - x3)
A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom
B = (x3^2 * (y1 - y2) + x2^2 * (y3 - y1) + x1^2 * (y2 - y3)) / denom
C = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom
В качестве альтернативы, если вы хотите выполнять математическую математику численно, вы, как правило, обращаетесь к системе линейной алгебры (например, ATLAS, хотя я не уверен, имеет ли она привязки C#/C++).
Спасибо Дэвид, я преобразовал ваш псевдокод в следующий код C#:
public static void CalcParabolaVertex(int x1, int y1, int x2, int y2, int x3, int y3, out double xv, out double yv)
{
double denom = (x1 - x2) * (x1 - x3) * (x2 - x3);
double A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom;
double B = (x3*x3 * (y1 - y2) + x2*x2 * (y3 - y1) + x1*x1 * (y2 - y3)) / denom;
double C = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom;
xv = -B / (2*A);
yv = C - B*B / (4*A);
}
Это то, что я хотел. Простой расчет вершины параболы. Я обработаю целочисленное переполнение позже.
Вы получаете следующие три уравнения путем прямой замены:
A*x1^2+B*x1+C=y1
A*x2^2+B*x2+C=y2
A*x3^2+B*x3+C=y3
Вы можете решить это, заметив, что это эквивалентно матричному произведению:
[x1^2 x1 1] [A] [y1]
|x2^2 x2 1|*|B| = |y2|
[x3^2 x3 1] [C] [y3]
Таким образом, вы можете получить A,B и C, инвертировав матрицу и умножив обратное на вектор справа.
Я вижу, что пока я публиковал эту статью, Джон Раш связал ее с учебным пособием, в котором более подробно рассматривается решение матричного уравнения, поэтому вы можете следовать этим инструкциям, чтобы получить ответ. Инвертировать матрицу 3х3 довольно легко, так что это не должно быть слишком сложно.
Я сделал нечто похожее на ответ @piSHOCK, также основанный на коде @AZDean. Если вам нужно запустить его интенсивно (или использовать его в Matlab, как я), это может быть самым быстрым.
Я предполагаю, что x1 == -1, x2 == 0, x3 == 1
,
a = y2 - ( y1 + y3) / 2 % opposite signal compared to the original definition of A
b = (y3 - y1) / 4 % half of the originally defined B
xExtr = b / a
yExtr = y2 + b * yExtr % which is equal to y2 + b*b / a
Вот код на Фортране, который реализует решение @david-z и @AZDean:
subroutine parabola_vertex(x1, y1, x2, y2, x3, y3, xv, yv)
real(dp), intent(in) :: x1, y1, x2, y2, x3, y3
real(dp), intent(out) :: xv, yv
real(dp) :: denom, A, B, C
denom = (x1 - x2) * (x1 - x3) * (x2 - x3)
A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom
B = (x3**2 * (y1 - y2) + x2**2 * (y3 - y1) + x1**2 * (y2 - y3)) / denom
C = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + &
x1 * x2 * (x1 - x2) * y3) / denom
xv = -B / (2*A)
yv = C - B**2 / (4*A)
end subroutine
def vertex(x1,x2,x3,y1,y2,y3):
'''Given three pairs of (x,y) points return the vertex of the
parabola passing through the points. Vectorized and common expression reduced.'''
#Define a sequence of sub expressions to reduce redundant flops
x0 = 1/x2
x4 = x1 - x2
x5 = 1/x4
x6 = x1**2
x7 = 1/x6
x8 = x2**2
x9 = -x7*x8 + 1
x10 = x0*x1*x5*x9
x11 = 1/x1
x12 = x3**2
x13 = x11*x12
x14 = 1/(x0*x13 - x0*x3 - x11*x3 + 1)
x15 = x14*y3
x16 = x10*x15
x17 = x0*x5
x18 = -x13 + x3
x19 = y2*(x1*x17 + x14*x18*x6*x9/(x4**2*x8))
x20 = x2*x5
x21 = x11*x20
x22 = x14*(-x12*x7 + x18*x21)
x23 = y1*(-x10*x22 - x21)
x24 = x16/2 - x19/2 - x23/2
x25 = -x17*x9 + x7
x26 = x0*x1*x14*x18*x5
x27 = 1/(-x15*x25 + y1*(x20*x7 - x22*x25 + x7) + y2*(-x17 + x25*x26))
x28 = x24*x27
return x28,x15 + x22*y1 + x24**2*x27 - x26*y2 + x28*(-x16 + x19 + x23)
Это пахнет домашней работой. "Спроси ученого" - это правильно. Скажите, что ваши 3 точки (x1, y1), (x2, y2) и (x3, y3). Затем вы получите три линейных уравнения:
| М11 М12 М13 | | A | | Z1 | | М21 М22 М23 | * | Б | = | Z2 | | М31 М32 М33 | | C | | Z3 |
Где M11 = x12, M12 = x1, M13 = 1, Z1 = y1 и аналогично для двух других строк, использующих (x2, y2) и (x3, y3) вместо (x1, у1).
Решение этой системы из 3 уравнений даст вам решение для A, B и C.
Работает на https://ideone.com/y0SxKU
#include <iostream>
using namespace std;
// calculate the vertex of a parabola given three points
// https://stackru.com/q/717762/16582
// @AZDean implementation with given x values
void CalcParabolaVertex(int x1, int y1, int x2, int y2, int x3, int y3, double& xv, double& yv)
{
double denom = (x1 - x2) * (x1 - x3) * (x2 - x3);
double A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom;
double B = (x3*x3 * (y1 - y2) + x2*x2 * (y3 - y1) + x1*x1 * (y2 - y3)) / denom;
double C = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom;
xv = -B / (2*A);
yv = C - B*B / (4*A);
}
// @piSHOCK immplementation assuming regular x values ( wrong!!! )
void CalcParabolaVertex2( int y1, int y2, int y3, double& xv, double& yv)
{
double d1 = y1 - y2;
double d2 = y1 - y3;
double a = -d1 + 0.5 * d2;
double b = 2 * d1 - 0.5 * d2;
double c = -y1;
xv = -0.5 * b / a;
yv = c - 0.25 * b * b / a;
}
// corrected immplementation assuming regular x values
void CalcParabolaVertex3( int y1, int y2, int y3, double& xv, double& yv)
{
double d1 = y1 - y2;
double d2 = y1 - y3;
double a = d1 - 0.5 * d2;
double b = -2 * d1 + 0.5 * d2;
double c = y1;
xv = -0.5 * b / a;
yv = c - 0.25 * b * b / a;
}
int main() {
double xv, yv;
CalcParabolaVertex( 0, 100, 1, 500, 2, 200, xv, yv );
cout << xv <<" "<< yv << "\n";
CalcParabolaVertex2( 100, 500, 200, xv, yv );
cout << xv <<" "<< yv << "\n";
CalcParabolaVertex3( 100, 500, 200, xv, yv );
cout << xv <<" "<< yv << "\n";
return 0;
}
Я добавил пару юнит-тестов для определения отрицательных пиковых значений: запуск в прямом эфире на https://ideone.com/WGK90S
Если вы сделаете предположение, что x1 == 0, x2 == 1, x3 == 2
Например, при локальной подгонке кривой к некоторому равномерно дискретизированному сигналу код @AZDean упрощается до:
d1 = y1 - y2
d2 = y1 - y3
a = d1 - 0.5 * d2
b = -2 * d1 + 0.5 * d2
c = y1
xVertex = -0.5 * b / a
yVertex = c - 0.25 * b * b / a