# Curve through four points -- simple algebra??

The motivation for this is Bezier curves. But, if you don't know what these are, you can skip down to the last paragraph, where the problem is described in purely algebraic terms.

Suppose I want to construct a quadratic Bezier curve that passes through some given points.

First, the standard approach: I would use three points, say $\mathbf A$, $\mathbf B$, $\mathbf C$, choose a "middle" parameter value, typically $h = \tfrac12$, and then construct a quadratic curve $\mathbf P(t)$, such that $\mathbf P(0) = \mathbf A$, $\mathbf P(h) = \mathbf B$, and $\mathbf P(1) = \mathbf C$. Easy.

But, it's fairly well known that given four coplanar points, you can (usually) construct a parabola that passes through them. This problem goes back to Newton, even. But, every parabola is a quadratic Bezier curve. So, given the four coplanar points $\mathbf A$, $\mathbf B$, $\mathbf C$, $\mathbf D$, I'm hoping that there is a fairly simple calculation that gives us a quadratic Bezier curve passing through them.

More specifically, let's call the quadratic Bezier curve $\mathbf P(t)$, again, and let $\mathbf P_0$, $\mathbf P_a$, $\mathbf P_1$ be its control points, so that $$\mathbf P(t) = (1-t)^2 \mathbf P_0 + 2t(1-t) \mathbf P_a + t^2 \mathbf P_1$$ We want to find control points $\mathbf P_0$, $\mathbf P_a$, $\mathbf P_1$ and parameter values $t=u$ and $t=v$ such that $$\mathbf P(0) = \mathbf A \quad ; \quad \mathbf P(u) = \mathbf B \quad ; \quad \mathbf P(v) = \mathbf C \quad ; \quad \mathbf P(1) = \mathbf D \quad$$ It's obvious that $\mathbf P_0 = \mathbf P(0) = \mathbf A$, and $\mathbf P_1 = \mathbf P(1) = \mathbf D$, so the remaining problem is to find the middle control point $\mathbf P_a$ and the two parameter values $u$ and $v$. I'd like to know how to calculate these remaining unknowns. Ideally a nice simple way.

There is a paper on this subject, but the authors' approach seems oddly complex and esoteric, to me. I'm hoping for a simple series of elementary calculations that a smart high-school student could understand.

As far as I can tell, you don't need to know anything about Bezier curves to solve this problem. It can be considered purely algebraically: we are given $\mathbf A$, $\mathbf B$, $\mathbf C$, $\mathbf D$ in $\mathbb R^2$, and we want to find $\mathbf P \in \mathbb R^2$ and $u,v \in \mathbb R$ such that $$(1-u)^2 \mathbf A + 2u(1-u) \mathbf P + u^2 \mathbf D = \mathbf B$$ $$(1-v)^2 \mathbf A + 2v(1-v) \mathbf P + v^2 \mathbf D = \mathbf C$$ So, it's really just an equation solving problem.

• I see now what you want to do! I have deleted my response. Commented May 20, 2013 at 10:21
• @Cantlog -- I was waiting for a better one. But I guess I've waited long enough :-) Commented Sep 27, 2013 at 0:00

Solved it myself, as below. I don't feel that this solution is "right" somehow, because it treats the interior points $$\mathbf B$$ and $$\mathbf C$$ very asymmetrically, so improvements are invited. But, anyway ...

Let's use the notation from the last few paragraphs of the question.

It's clear that $$\mathbf P(0) = \mathbf A$$, and $$\mathbf P(1) = \mathbf D$$, so these two points are interpolated, already, and we only have to worry about the other two points, $$\mathbf B$$ and $$\mathbf C$$.

First, we find numbers $$h$$ and $$k$$ such that $$\mathbf C = \mathbf A + h(\mathbf B - \mathbf A) + k(\mathbf D - \mathbf A) = (1-h-k)\mathbf A + h\mathbf B + k\mathbf D$$ This is possible provided that $$\mathbf A$$, $$\mathbf B$$, $$\mathbf D$$ are not collinear. Then, since $$\mathbf P(v) = \mathbf C$$, we have $$(1-h-k)\mathbf A + h\mathbf B + k\mathbf D = (1-v)^2 \mathbf A + 2v(1-v) \mathbf P + v^2 \mathbf D$$ But, since $$\mathbf P(u) = \mathbf B$$, we know that $$\mathbf B = (1-u)^2 \mathbf A + 2u(1-u) \mathbf P + u^2 \mathbf D$$ Substituting for $$\mathbf B$$ on the left-hand side, and equating coefficients of $$\mathbf A$$, $$\mathbf P$$, $$\mathbf D$$ gives $$(1-v)^2 = 1 - h - k +h(1-u)^2$$ $$2v(1-v) = 2hu(1-u)$$ $$v^2 = hu^2 + k$$ We can easily eliminate $$v$$ from these last three equations using the fact that $$[2v(1-v)]^2 = 4[v^2][(1-v)^2]$$. We get: $$[2hu(1-u)]^2 = 4[hu^2 + k][1 - h - k +h(1-u)^2]$$ After a little algebra, this reduces to: $$h(1-h)u^2 - 2hku + k(1-k) = 0$$

So, we solve this quadratic for $$u$$, and then get the unknown interior control point $$\mathbf P$$ from

$$\mathbf P = \frac{\mathbf B - (1-u)^2 \mathbf A -u^2 \mathbf D}{2u(1-u)}$$

The number of solutions depends on the number of real solutions of the quadratic. Quite often, you can draw two different quadratic Bezier curves through the four given points, as explained in the paper cited in the question.

The quoted paper is actually quite smart, not complex or esoteric.

I wrote some C++ code, using Eigen, to find the parabola's if they exist, using their method. Here Q1 through Q4 are the four points whose coordinates can be accessed with the member functions x() and y().

      Eigen::Matrix3d R;
R << Q1.x(), Q1.y(), 1.0,
Q2.x(), Q2.y(), 1.0,
Q3.x(), Q3.y(), 1.0;

Eigen::RowVector3d Q4_1;
Q4_1 << Q4.x(), Q4.y(), 1.0;

auto q = Q4_1 * R.inverse();


Now we have the vector q that can be used to see if there is a solution:

      int count = 0;
for (int i = 0; i < 3; ++i)
{
// No solution if two q_i values are equal to 1.
if (utils::almost_equal(q(i), 1.0, 10e-9))
++count;
}
if (count >= 2 || q(0) * q(1) * q(2) >= 0.0)
{
// There is no parabola going through all four points.
continue;
}


If we get here there are one or two solutions. I didn't bother to check for that (but well - one if the discriminant is zero).

      for (int solution = 0; solution < 2; ++solution)
{
double discriminant = utils::square(-2.0 * q(1) * q(2)) - 4.0 * q(1) * (1.0 - q(1)) * q(2) * (1.0 - q(2));
// Calculate alpha.
double a = (2.0 * q(1) * q(2) + ((solution == 0) ? -1 : 1) * std::sqrt(discriminant)) / (2.0 * q(1) * (1.0 - q(1)));

Eigen::Matrix3d V;
V << 0.0, 0.0, 1.0,
a*a,   a, 1.0,
1.0, 1.0, 1.0;

Eigen::Matrix3d V_inverse = V.inverse();


At this point you can plot the parabola defined by V_inverse * R. For example,

        std::vector<Point> curve_points;
{
for (int i = -200; i <= 200; ++i)
{
double t = i * 0.1;
Eigen::RowVector3d T;
T << (t * t), t, 1.0;
auto Phi = T * (V_inverse * R);
curve_points.emplace_back(Phi(0), Phi(1));
}
}
// Next solution.
}