[Math] Cylinder in 3D from five points

geometrylinear algebranonlinear system

I have five points in 3D space which I want to use to form a cylinder in 3D. I assume that these points are exactly on the cylinder.

I know the equation for a cylinder is given by an origin and a normal to define the center axis of the cylinder, along with the cylinder's radius.

The solution I came up with was to project the points onto a 2D plane via two unnormalized vectors $u_{ox}$ and $u_{oy}$ that are orthogonal to each other and to the normal, and then fitting the points to a circle equation.

$(p_i \cdot u_{0x} – c_x)^2 + (p_i \cdot u_{0y} – c_y)^2 = r^2$

Depending on your choice of $u_{0x}$ and $u_{0y}$, when you expand the equations you get the solution (for paticular $a,b,…,h$)

$x^2 + a \cdot x \cdot y + b \cdot x \cdot z + c \cdot x + d \cdot y^2 + e \cdot y \cdot z + f \cdot y + g \cdot z^2 + h = 0$

(Is this called a generic conic equation?)

Now I know it would take nine sets of points $[x,y,z]=p_i$ to find the solution to the nine coefficients, but I have additional knowledge of the relationship between the coefficients.

Going back to the original circle equation, you can choose some nice unnormalized vectors to develop your coefficients (provided the normal is not simply $+x$).

$u_{0x}=\left[\begin{array}{c c c}
1 & -u_x \cdot u_y & u_x
\end{array}\right]$

$u_{0y}=\left[\begin{array}{c c c}
0 & s & u_y \cdot s
\end{array}\right];
s = \sqrt{\frac{1+u_x^2\cdot u_y^2+u_x^2}{1+u_y^2}}$

While $u_{ox}$ and $u_{oy}$ are orthogonal to each other, it requires a scaling factor $s$ or you'll end up with an ellipse instead of a circle when you do your projection.

Between the circle equation and the equation for the unnormalized vectors, I've got five unknowns $(u_x, u_y, c_x, c_y, r)$ and the scaling factor $s$ (for which I have a gnarly equation). That should mean that only five sets of points $p_i$ are necessary to find what these unknowns are. But if you only work on the generic conic equation, you're short by four equations to find all of the unknowns.

In terms of what all of those coefficients look like…

$\begin{array}{r\;l}
p_x^2 & ( 1 ) \\
+ p_x \cdot py & ( -2 u_x \cdot u_y) \\
+ p_x \cdot pz & ( -2 u_x ) \\
+ p_x & ( -2 cx) \\
+ p_y^2 & (u_x^2 \cdot u_y^2 + s^2) \\
+ p_y \cdot p_z & ( -2 u_x^2 u_y +2 s^2 u_y ) \\
+ p_y & ( 2 u_x \cdot u_y \cdot c_x – 2 s \cdot c_y) \\
+ p_z^2 & ( u_x^2 + s^2 \cdot u_y^2)\\
+ p_z & ( 2 u_x \cdot c_x – 2 s \cdot u_y \cdot c_y) \\
+ & (c_x^2 + c_y^2 – r^2) \\
= & 0
\end{array}$

If I knew $a, b, …, h$, I could easily find $u_x, u_y$, etc., and then translate them into a more proper $origin$, $normal$, and $radius$, but I'm having a math brain fart. How do I get $a, b, …, h$ with only five points?

Anyone got some thoughts on how to proceed? Or know of a simpler solution?

Best Answer

Interesting question!

(Is this called a generic conic equation?)

To me, a generic conic would be something planar. You essentially describe the 3d analogon here, so I'd call this a quadric.

I know it would take nine sets of points $[x,y,z]=p_i$ to find the solution to the nine coefficients

And yet you only use 8 variables $a\ldots h$. Well, as you found out and wrote in a comment, that's because you forgot the linear term in $z$. So yes, 9 points would match 9 variables. You assume (without explicitely stating that assumption) a constant coefficient of $1$ for the $x^2$ term. That's permissible except if that coefficient would turn out to be zero. Even more general, you'd take a tenth variable for the $x^2$ term, then use 9 points to determine said variables up to a common scalar factor which is irrelevant as the right hand side is zero. The equation is homogeneous.

Anyone got some thoughts on how to proceed? Or know of a simpler solution?

Abbreviating $x_i:=p_i \cdot u_{0x}$ and $y_i:=p_i \cdot u_{0y}$ as the projected in-plane coordinates of the point, I'd write the circle as

$$(x_i^2+y_i^2)+ax_i+by_i+c=0$$

That way you have 6 unknowns: $u_x,u_y,s,a,b,c$. This description feels easier than your parametrization using $c_x,c_y,r$. But the relations between these are still non-linear, which makes the whole problem quite tricky. You probably should aim to make all conditions polynomial, i.e. rewrite your condition for $s$ to

$$s^2(1+u_y^2)=1+u_x^2\cdot u_y^2+u_x^2$$

You can of course try to let some computer algebra system solve this for you. Or you can look in to elimination theory, and start eliminating some of the unknowns from your equations using Gröbner bases or resultants. Which is probably what your computer algebra system would do internally.

I've actually tried to do the CAS computation with a bunch of random points. Wolfram Alpha doesn't accept as much input as I'd need for this, so I'll be using Sage, and in particular its method for computing the variety of an ideal.

At first I had made a mistake when I tried to use your $u_{0x}$ and $u_{0y}$, see later section for details. So looking for an alternate solution I defined $u_{0x}=(x_1,y_1,0)$ and $u_{0y}=(x_2,y_2,z_2)$ and then added conditions for $u_{0x}\cdot u_{0x}=1$, $u_{0y}\cdot u_{0y}=1$ and $u_{0x}\cdot u_{0y}=0$. Or in other words, both vectors should be unit length and orthogonal to one another. I've assumed that the $z$ coordinate of the axis is non-zero, while you chose $x\neq 0$, but that shouldn't matter. (If you automate this computation, handling such special cases automatically can be troublesome.)

A specific example

With this formulation I tried to solve a specific instance in Sage. I used five points, chosen at random from $\{0,1,2,3,4\}^3$.

$$ p_1=(1, 4, 4)\quad p_2=(4, 3, 4)\quad p_3=(2, 0, 1)\quad p_4=(0, 0, 4)\quad p_5=(2, 4, 0) $$

For these I got $24$ distinct solutions. Actually I usually got $24$ solutions over $\mathbb C$ but only some multiple of $8$ of them real, so the above example was specifically chosen from a number of random examples as one where all the solutions are real.

Here is the Sage code to compute these solutions:

# pts = random_matrix(ZZ, 5, 3, x=0, y=5).rows() # for new random instance
pts = matrix([(1,4,4), (4,3,4), (2,0,1), (0,0,4), (2,4,0)]).rows() # use my case
PR.<a,b,c,x1,y1,z1,x2,y2,z2> = QQ[] # multi-variate polynomial ring
u1 = vector((x1,y1,z1))
u2 = vector((x2,y2,z2))
id1 = ideal([(p*u1)^2 + (p*u2)^2 + a*p*u1 + b*p*u2 + c for p in pts]
          + [z1-0, u1*u2, u1*u1-1, u2*u2-1])
id1.dimension() # must be zero
len(id1.variety(CC)) # will be 24 for most examples
len(id1.variety(RR)) # will be 24 for my example, and usually a multiple of 8

These solutions have to occur in groups of four, since one can negate all coordinates of $u_{0x}$ and/or $u_{0y}$ and still describe essentially the same conic. But this still leaves $24/4=6$ really different cylinders for the situation above. Their axis directions according to my computation using the example above are approximately

$$ (0.678, 0.729, 0.0989) \quad (-0.605, 0.489, -0.628) \quad (-0.408, 0.179, -0.895) \\ (-0.422, -0.104, 0.900) \quad (-0.0774, -0.964, 0.253) \quad (-0.963, -0.108, 0.245) $$

These six coordinates also correspond to six different values, e.g. for the unknown $c$. These six values are the six roots of the polynomial

\begin{align*} 16947006302936317345&\cdot c^6\\ - 423964877845474796312&\cdot c^5\\ + 3414707935712257584272&\cdot c^4\\ - 8163333268350816357120&\cdot c^3\\ - 9034354158004982553600&\cdot c^2\\ + 35770219295155494912000&\cdot c\\ + 13673321316256573440000&= 0 \end{align*}

This I found by exact computation using algebraic numbers (QQbar in Sage, but note that the computation can take considerable time). The Galois group of this polynomial is $S_6$, which is not solvable. So according to Galois theory there can be no way to express the value of $c$ for any of these solutions using radicals. Which in turn means that I'd strongly advise against solving this kind of problem without help from a computer to do the heavy number crunching.

Back to your set of variables

When I first used your vhoice of $u_{0x}$ and $u_{0y}$, I must have forgotten some extra variables, which made me believe that there were additional real degrees of freedom which I couldn't solve for. Once I realized and fixed that mistake, I could apply the above approach to your parametrization as well. I found $12$ solutions, instead of my $24$, which is to be expected as you fixed one coordinate to $1$ which fixed the sign of that vector. So my choice of $u_{0x}$ and $u_{0y}$ is probably not better than yours, while the parametrization of the circle using just three linear unknown variables should still be a good idea.

Feeding my example into your description, the variables that define your vectors are characterized by these polynomials:

\begin{align*} \small 1478610113817600&\cdot u_x^6 & \small 266500&\cdot u_y^6 & \small 2186287868683696026446069760000&\cdot s^{12} \\ \small - 3863371808155392&\cdot u_x^5 & \small + 539357&\cdot u_y^5 & \small -33034531469518995706122669195264&\cdot s^{10} \\ \small - 3541905515941632&\cdot u_x^4 & \small - 8742480&\cdot u_y^4 & \small +96487595035479832160430206853120&\cdot s^8 \\ \small + 565191711017856&\cdot u_x^3 & \small + 12151234&\cdot u_y^3 & \small -106608759197577962918341435697664&\cdot s^6 \\ \small + 691347653169852&\cdot u_x^2 & \small - 5732364&\cdot u_y^2 & \small +44504844774915504691341808859664&\cdot s^4 \\ \small + 59575310863341&\cdot u_x & \small + 1024565&\cdot u_y & \small -3621892801242489875603841840033&\cdot s^2 \\ \small - 1445045015375 &= 0 & \small - 59000 &= 0 & \small + 75225907399919929034745390625 &= 0 \end{align*}

So you “only” need to compute the roots and check which combinations fit together, then you know how to project the points and continue from there. Perhaps this is of use at some point, if only to illustrate how hard an exact solution can be to describe even for very simple input.

Related Question