Cylinder Fitting – Elite Math Program

Several years ago, I was given a practical geometry problem:

Given thousands of measured points in three-dimensional space, find the cylinder that best fits them.

There was an additional requirement that made the problem much more interesting. The calculation had to be fast enough to operate in real time.

This meant that simply defining a mathematically reasonable fitting problem was not enough. The equations had to be arranged so that thousands of data points could be processed efficiently, the derivatives could be evaluated quickly, and each optimization step would require solving only a very small system of equations.

The resulting method combines geometry, calculus, least squares, and the Gauss–Newton method.

The Geometry of a Cylinder

Suppose that

X = {(x_{1}, x_{2}, x_{3})}^{T}

is a point on the axis of a cylinder, and

U = {(u_{1}, u_{2}, u_{3})}^{T}

is a vector pointing along the axis.

Let $P$ be a point on the surface of the cylinder.

The defining geometric property is simple:

The perpendicular distance from every point on the cylinder to its axis is the radius $r$ .

Distance from a Point to the Axis

Consider the vector

P - X .

Part of this vector points along the cylinder axis, while the remaining part is perpendicular to the axis.

If $α$ is the angle between $P - X$ and $U$ , then

‖ P - X ‖ sin α = r .

Squaring gives

{‖ P - X ‖}^{2} {sin α}^{2} = r^{2} .

The dot product gives

{({(P - X)}^{T} U)}^{2} = {‖ P - X ‖}^{2} {‖ U ‖}^{2} {cos α}^{2} .

Combining these relations eliminates the angle completely and gives

{‖ P - X ‖}^{2} - \frac{{({(P - X)}^{T} U)}^{2}}{{‖ U ‖}^{2}} = r^{2} .

This is the equation we need.

The first term measures the squared distance from $X$ to $P$ . The fraction subtracts the squared component parallel to the axis. What remains is exactly the squared perpendicular distance to the axis.

Turning the Geometry into a Residual

Measured points will not lie exactly on a perfect cylinder.

For a measured point $P$ , define

R (P) = {‖ P - X ‖}^{2} - \frac{{({(P - X)}^{T} U)}^{2}}{{‖ U ‖}^{2}} - r^{2} .

For a point exactly on the cylinder,

R (P) = 0 .

For a noisy measured point, the residual will usually be nonzero.

Reducing the Number of Unknowns

At first the cylinder appears to have seven parameters: three coordinates for a point on its axis, three coordinates for the axis direction, and the radius.

But there is redundancy.

Moving the reference point $X$ along the axis does not change the cylinder. Also, multiplying $U$ by a nonzero constant does not change its direction.

Assuming the cylinder axis is not parallel to the $(x_{1}, x_{2})$ plane, choose

X = {(x_{1}, x_{2}, 0)}^{T}, U = {(u_{1}, u_{2}, 1)}^{T} .

Now the entire cylinder is described by only five unknown numbers:

S = {(x_{1}, x_{2}, u_{1}, u_{2}, r)}^{T} .

Two numbers locate the axis, two determine its direction, and one gives the radius.

From Geometry to Least Squares

Suppose our measuring system produces $m$ points:

P_{1}, P_{2}, \dots, P_{m} .

Each point produces a residual $R_{j}$ .

We find the best-fitting cylinder by minimizing

f (S) = \frac{1}{2} \sum_{j = 1}^{m} R_{j}^{2} .

This is a nonlinear least-squares problem.

Why the Derivatives Matter

Because the calculation had to run quickly, I derived the residual derivatives analytically.

For each measured point, we need derivatives with respect to the five unknown parameters:

\nabla R = (\frac{\partial R}{\partial x_{1}}, \frac{\partial R}{\partial x_{2}}, \frac{\partial R}{\partial u_{1}}, \frac{\partial R}{\partial u_{2}}, \frac{\partial R}{\partial r}) .

The formulas contain several repeated expressions. In an implementation, there is no reason to calculate the same quantity over and over.

For example, define

q = \frac{(p_{1} - x_{1}) u_{1} + (p_{2} - x_{2}) u_{2} + p_{3}}{u_{1}^{2} + u_{2}^{2} + 1} .

Then several derivative formulas become much shorter and faster to evaluate.

This may look like a minor algebraic simplification on paper. In a real-time numerical algorithm that evaluates the same expressions thousands of times, such simplifications matter.

The Jacobian

Collect all the residuals into the vector

E = {(R_{1}, R_{2}, \dots, R_{m})}^{T} .

The Jacobian $J$ has one row for every measured point and one column for every unknown parameter.

Therefore,

J is an m \times 5 matrix.

Gauss–Newton

The least-squares objective can be written compactly as

f = \frac{1}{2} {‖ E ‖}^{2} .

Its gradient is

\nabla f = J^{T} E .

Gauss–Newton approximates the Hessian by

\nabla^{2} f \approx J^{T} J .

At each iteration, instead of solving the original nonlinear problem from scratch, we solve

J^{T} J G = - J^{T} E .

The vector $G$ tells us how to change the current estimate of the cylinder.

Why Thousands of Points Are Not as Bad as They Sound

This is the computational feature that makes the method especially attractive.

Suppose the measuring system gives us 10,000 points.

Then $J$ has 10,000 rows.

But it still has only five columns.

Therefore,

J^{T} J

is only a

5 \times 5

matrix.

And

J^{T} E

contains only five numbers.

So although every iteration uses information from thousands of measured points, the linear system that determines the next step has only five unknowns.

That is a very useful structure for a real-time calculation.

Do Not Always Take the Full Step

A Gauss–Newton direction tells us which way to move, but taking the entire step is not always wise.

Write the update as

S_{k + 1} = S_{k} + a_{k} G_{k}, 0 < a_{k} \leq 1 .

The step length is chosen adaptively. Start with $a_{k} = 1$ and reduce it if necessary until the new point produces a sufficient decrease in the objective function.

In practice, the full step can often be accepted, but the line search gives the algorithm additional protection when the current estimate is not yet close to the solution.

When Do We Stop?

At a minimum we expect

\nabla f = 0 .

Numerically, we stop when

‖ \nabla f ‖ < ε,

where $ε$ is a small tolerance.

There is no benefit in demanding far more numerical precision than the data and the computer arithmetic can support. An unnecessarily small tolerance can even create numerical difficulties.

Did It Work?

Yes.

This was not a numerical example invented after the fact. The cylinder fitting problem came from a real application, and the algorithm had to work fast enough for real-time use.

The implementation worked extremely well.

What I found especially satisfying was that the final algorithm was built from familiar mathematical ingredients:

geometry → least squares → calculus → Gauss–Newton → real-time computation.

The important part was arranging those ingredients in the right way.

From Cylinders to Cones

A cylinder has constant radius. Once this fitting method works, a natural question is:

What happens if the radius changes as we move along the axis?

That leads to a cone.

The cone-fitting problem requires one additional parameter—the opening angle—and the derivatives become more complicated. But the central strategy remains the same:

choose an efficient geometric residual → derive its Jacobian → form a nonlinear least-squares problem → solve it with Gauss–Newton.

That will be the subject of the next post.

The Larger Lesson

There is an important distinction between solving a mathematical problem and solving it in a form that is useful in practice.

With 10,000 measured points, the original problem sounds large. But the geometry allows the unknown cylinder to be represented by only five parameters. Gauss–Newton then converts each iteration into a small five-variable linear problem.

This is one of the recurring ideas in applied mathematics:

A good mathematical formulation can be as important as the algorithm itself.

When the formulation is right, a problem involving thousands of three-dimensional measurements can become small enough to solve in real time.

A Closely Related Fitting Problem

Cylinder fitting is only one version of the real-time surface-fitting problem. A closely related challenge arises when the measured surface is a cone rather than a cylinder.

The change in geometry leads to a different mathematical problem, but the objective is the same: use the structure of the surface to obtain an accurate fit quickly enough for practical real-time computation.

Continue exploring: How Do You Fit a Cone to Thousands of Points in Real Time?

Tag: Cylinder Fitting

How Do You Fit a Cylinder to Thousands of Points in Real Time?