Elite Math Program

Tag: 3D Geometry

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

    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 = ( x1 , x2 , x3 ) T

    is a point on the axis of a cylinder, and

    U = ( u1 , u2 , u3 ) 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 = r2.

    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 − ( (P−X) T U ) 2 ‖U‖ 2 = r2 .

    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 − ( (P−X) T U ) 2 ‖U‖ 2 − r2 .

    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 ( x1 , x2 ) plane, choose

    X = ( x1 , x2 , 0 ) T , U = ( u1 , u2 , 1 ) T .

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

    S = ( x1 , x2 , u1 , u2 , 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:

    P1 , P2 , … , Pm .

    Each point produces a residual Rj .

    We find the best-fitting cylinder by minimizing

    f(S) = 12 ∑ 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:

    ∇R = ( ∂R ∂x1 , ∂R ∂x2 , ∂R ∂u1 , ∂R ∂u2 , ∂R ∂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 = (p1 −x1) u1 + (p2 −x2) u2 + p3 u12 + u22 + 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 = ( R1 , R2 , … , Rm ) T .

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

    Therefore,

    J is an m×5 matrix.

    Gauss–Newton

    The least-squares objective can be written compactly as

    f = 12 ‖E‖ 2 .

    Its gradient is

    ∇f = JTE .

    Gauss–Newton approximates the Hessian by

    ∇2f ≈ JTJ .

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

    JT J G = − JT 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,

    JTJ

    is only a

    5×5

    matrix.

    And

    JTE

    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

    Sk+1 = Sk + ak Gk , 0 < ak ≤ 1 .

    The step length is chosen adaptively. Start with ak=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

    ∇f = 0 .

    Numerically, we stop when

    ‖∇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?

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

    Several years ago, I was given an applied mathematics problem that sounded simple at first:

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

    But there was an additional requirement that changed the problem completely. The calculation had to be fast enough to be used in real time.

    It is one thing to describe an optimization problem mathematically. It is another thing to solve it repeatedly, quickly, and reliably while a real system is operating.

    The solution came from combining elementary geometry with nonlinear least squares and the Gauss–Newton method.

    The Geometric Problem

    Suppose a collection of measured points

    P1 , P2 , … , Pm

    lies close to the surface of an unknown right circular cone.

    We want to determine three things:

    the vertex of the cone, the direction of its axis, and its opening angle.

    Let

    X = ( x1 , x2 , x3 )

    be the vertex, let U point along the axis, and let α be the half-angle of the cone.

    If P lies exactly on the cone, then the vector

    P−X

    makes angle α with the cone axis.

    The dot-product formula for the angle between two vectors therefore gives

    ( (P−X) T U ) 2 = ‖P−X‖ 2 ‖U‖ 2 cosα2 .

    This one equation contains the geometry of the cone.

    Turning Geometry into an Error

    Measured data will not lie exactly on a perfect cone. There will be noise, measurement error, and small deviations from the ideal surface.

    So instead of asking whether a point satisfies the cone equation exactly, we measure how far it is from satisfying the equation.

    Define the residual

    R (P) = ‖P−X‖ 2 cosα 2 − ( (P−X) T U ) 2 ‖U‖ 2 .

    For a point exactly on the cone,

    R(P) =0.

    For a measured point near the cone, the residual will generally be nonzero.

    Why Use This Residual?

    There is an important practical point here.

    One could try to calculate the exact shortest geometric distance from every measured point to the cone. But when thousands of points must be processed repeatedly in real time, the computational form of the problem matters.

    The residual above is obtained directly from the cone equation. It can be evaluated using additions, multiplications, dot products, and a few simple functions.

    This was exactly what was needed for the application: a mathematical description that could be differentiated explicitly and evaluated rapidly.

    Six Unknown Parameters

    The direction of the axis does not depend on the length of U. Multiplying the axis vector by a nonzero constant gives the same axis.

    Assuming the axis is not parallel to the (x1, x2) plane, we can therefore write

    U = ( u1 , u2 , 1 ) .

    The entire cone is then described by only six unknown numbers:

    S = ( x1 , x2 , x3 , u1 , u2 , α ) T .

    Three numbers locate the vertex, two determine the axis direction, and one determines the cone angle.

    From One Point to Thousands of Points

    Suppose there are m measured points.

    For each point Pj , compute a residual Rj .

    We then look for the cone that minimizes the total squared residual:

    f(S) = 12 ∑ j=1 m R j 2 .

    This is a nonlinear least-squares problem.

    Why Not Just Solve the Equations?

    If the measurements were perfect, six carefully chosen points might appear to be enough to determine six unknown parameters.

    Real data do not work that way.

    Measurements contain noise, and a small collection of points may give a poor estimate. Instead, we can use hundreds or thousands of points simultaneously and find the cone that best fits all of them.

    But now the equations are nonlinear. There is no simple matrix formula that immediately gives the answer.

    This is where Gauss–Newton enters the story.

    The Jacobian

    Put all the residuals into one vector:

    E = ( R1 , R2 , … , Rm ) T .

    Then

    f = 12 ‖E‖ 2 .

    The Jacobian J contains the derivatives of every residual with respect to the six cone parameters.

    Its j-th row is

    Jj = ( ∂Rj ∂x1 , ∂Rj ∂x2 , ∂Rj ∂x3 , ∂Rj ∂u1 , ∂Rj ∂u2 , ∂Rj ∂α ) .

    In the actual implementation, these derivatives were derived analytically rather than estimated numerically.

    That requires more work at the beginning, but once the formulas are known, they can be evaluated very efficiently.

    The Key Gauss–Newton Approximation

    The gradient of the least-squares objective has a particularly simple form:

    ∇f = JT E .

    The exact Hessian is

    ∇2f = JTJ + ∑ j=1 m Rj ∇2 Rj .

    Computing the second term repeatedly is considerably more expensive.

    Gauss–Newton makes the approximation

    ∇2f ≈ JTJ .

    This approximation becomes especially natural when the fit is already good, because the residuals Rj are then small.

    That observation was crucial for a real-time implementation. Instead of constructing the complete second derivative at every iteration, we can work primarily with the Jacobian.

    One Iteration

    At the current estimate Sk , the Gauss–Newton direction Gk is obtained by solving

    JkT Jk Gk = − JkT Ek .

    Then update the six cone parameters:

    Sk+1 = Sk + ak Gk .

    Here 0<ak≤1 is the step length.

    Why Use a Line Search?

    Taking the full Gauss–Newton step is not always a good idea when the current estimate is still far from the solution.

    The implementation therefore reduces the step when necessary until the objective function decreases sufficiently.

    In mathematical form, choose the step length so that

    f ( Sk + ak Gk ) ≤ f (Sk) + ρ ak ∇fk T Gk .

    The purpose is simple: do not accept an iteration that moves too aggressively in a direction that fails to improve the fit sufficiently.

    The Real-Time Idea

    At first glance this may still look like a large calculation. There may be thousands of measured points.

    But notice something important.

    No matter how many data points we have, there are still only six unknown cone parameters.

    The Jacobian may have thousands of rows, but only six columns. The matrix

    JTJ

    is therefore only

    6×6.

    Similarly,

    JTE

    has only six components.

    So the data set can be large while the system that must be solved at each Gauss–Newton iteration remains very small.

    This is exactly the type of structure one wants to exploit in a real-time numerical algorithm.

    Did It Work?

    Yes.

    This was not an exercise invented to demonstrate Gauss–Newton. I developed the method because the cone-fitting calculation was needed in a real application, and it had to operate in real time.

    The implementation worked extremely well.

    That experience taught me an important lesson about applied mathematics: the best mathematical formulation is not necessarily the one that looks most sophisticated on paper.

    Sometimes the decisive question is:

    Can we formulate the problem so that the computer can solve it quickly enough to be useful?

    Why the Approximation Works

    There is another nice feature of Gauss–Newton.

    Recall that the exact Hessian is

    ∇2f = JTJ + ∑ j=1 m Rj ∇2 Rj .

    As the estimated cone approaches the data, the residuals become small. Consequently, the second term becomes less important and

    ∇2f ≈ JTJ .

    So as the algorithm approaches a good fit, the inexpensive approximation becomes increasingly appropriate.

    A Final Check

    Although computing the full Hessian at every iteration is expensive, it can still be useful after the optimization has finished.

    At the final solution we can evaluate

    ∇2f = JTJ + ∑ j=1 m Rj ∇2 Rj .

    If its eigenvalues are positive, the Hessian is positive definite and the computed stationary point is a local minimum.

    In other words, we use the inexpensive approximation while speed matters, and we can use the more expensive calculation afterward as a check.

    The Larger Lesson

    The mathematics of this problem can be summarized in one chain:

    3D measurements  →  cone geometry  →  residuals  →  least squares  →  Jacobian  →  Gauss–Newton  →  real-time fit.

    The cone itself is elementary geometry. Least squares is a classical idea. The derivatives require calculus. Gauss–Newton comes from numerical optimization.

    None of these ingredients alone solves the practical problem.

    The solution comes from putting them together in a form that a computer can evaluate rapidly.

    That is one of the most satisfying aspects of applied mathematics: sometimes a few familiar mathematical ideas, combined in the right way, turn a difficult real-world computation into something that works almost like magic.

    A Closely Related Fitting Problem

    Fitting a cone to a large cloud of measured points is one example of a broader problem: how can we recover a geometric surface accurately when the computation has to be fast enough for real-time use?

    A closely related problem replaces the cone by a cylinder. The geometry changes, but the same practical challenge remains: turn thousands of measurements into a reliable geometric fit without an expensive general-purpose optimization.

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