A non-iterative method for robustly computing the intersections between a line and a curve or surface

by   Xiao Xiao, et al.
University of Cambridge

The need to compute the intersections between a line and a high-order curve or surface arises in a large number of finite element applications. Such intersection problems are easy to formulate but hard to solve robustly. We introduce a non-iterative method for computing intersections by solving a matrix singular value decomposition (SVD) and an eigenvalue problem. That is, all intersection points and their parametric coordinates are determined in one-shot using only standard linear algebra techniques available in most software libraries. As a result, the introduced technique is far more robust than the widely used Newton-Raphson iteration or its variants. The maximum size of the considered matrices depends on the polynomial degree q of the shape functions and is 2q × 3q for curves and 6 q^2 × 8 q^2 for surfaces. The method has its origin in algebraic geometry and has here been considerably simplified with a view to widely used high-order finite elements. In addition, the method is derived from a purely linear algebra perspective without resorting to algebraic geometry terminology. A complete implementation is available from http://bitbucket.org/nitro-project/.



There are no comments yet.


page 1

page 2

page 3

page 4


On the number of intersection points of lines and circles in ℝ^3

We consider the following question: Given n lines and n circles in ℝ^3, ...

Bounds for Substituting Algebraic Functions into D-finite Functions

It is well known that the composition of a D-finite function with an alg...

Changing Views on Curves and Surfaces

Visual events in computer vision are studied from the perspective of alg...

An adaptive high-order surface finite element method for the self-consistent field theory on general curved surfaces

In this paper, we develop an adaptive high-order surface finite element ...

Determining surfaces of revolution from their implicit equations

Results of number of geometric operations (often used in technical pract...

Interrogation of spline surfaces with application to isogeometric design and analysis of lattice-skin structures

A novel surface interrogation technique is proposed to compute the inter...

Spectral mesh-free quadrature for planar regions bounded by rational parametric curves

This work presents spectral, mesh-free, Green's theorem-based numerical ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

There has recently been an increased academic and industrial interest in high-order finite elements due to their efficiency advantages over classical low-order elements, see e.g.[1, 2]. To achieve their full potential high-order methods require the curved domain boundaries to be approximated with non-planar elements. The intersection between the curved elements and a given line is required in a wide range of applications, including contact [3], mesh generation [4, 5] and immersed finite elements [6]. As known, intersection computations lead to an easy to formulate, but hard to solve, nonlinear root-finding problem. The prevalent technique in computational mechanics for solving such problems is Newton-Raphson iteration, which is in general not very robust.

In computer-aided geometric design (CAD) and manufacturing (CAM) intersection computation is a recurring task and, to this end, a number of ingenious methods have been developed [7]. Especially promising are the non-iterative methods with a origin in algebraic geometry, which are for the most part unknown in computational mechanics. Algebraic geometry deals with systems of polynomial equations and geometric objects defined by them and provides the most rigorous framework for intersection computations [8]. The specific technique considered in this paper is the implicisation technique proposed by Busé [9], which shares some commonalities with the method of moving lines/planes introduced by Sederberg and Chen [10]. Different from the original work, in this paper we derive the method from a purely linear algebra viewpoint. To follow the presented derivations it is sufficient to only know the notion of the algebraic degree of a curve or surface. The algebraic degree of a curve or surface is defined as its number of intersections with a line. Counting all intersections (real, complex, multiple ones and ones at infinity), the algebraic degree of a polynomial curve of degree 

and the corresponding tensor-product surface are 

and , respectively. As will become clear, adopting a linear algebra viewpoint has the added benefit that many linear algebra techniques, like preconditioning and various matrix decompositions, become readily available for intersection computations.

In the following, we first discuss the intersection between a planar Lagrange curve and a line and provide an easy to follow illustrative example. This simple case is sufficient to introduce and discuss the key aspects of the proposed non-iterative technique. Its extension to the surface case is straightforward and is discussed last.

2 Intersection of lines with curves

2.1 Moving lines and intersections

Let  be a planar parametric curve, with , of degree  given either in Lagrange basis  or power (monomial) basis  with


where  are the nodal coordinates and  are the coefficients in the power basis. As usual, the power basis  contains the consecutive powers of  from  up to . The two basis are related by


where  is the Vandermonde matrix and  is the parametric coordinate of the -th Lagrange node. Following a similar approach a curve given in any other polynomial basis can be re-expressed in the power basis.

To define a point  on the curve as the intersection of several moving lines, or pencils of lines, consider



is an auxiliary vector collecting the parameters of the line. For a fixed 

equation (3) describes a line and the line moves with the parameter , see Figure (a)a . The three parameters of the line are assumed to be polynomial functions given by


The degree  of the power basis  has to be chosen sufficiently high in order to be able to compute all the intersection points (real, complex, multiple ones and ones at infinity). The number of intersection points is equivalent to the algebraic degree of the curve . A curve  of degree  has  intersection points with a line. As will become clear, the number of intersections implies a constraint on the minimum possible degree  for .

(a) Two linear moving lines () defining a quadratic curve ()
(b) Four of the five cubic moving lines () for a cubic curve () at and
Figure 3: Moving lines and their intersections.

Next, the aim is to find several lines (3), or more specifically their coefficients  in (4), with a common intersection point which is a point on the curve . It is required that each line satisfies at the common intersection point


After introducing the definition of the curve (1) this yields


The bracketed term can be expressed in a new power basis  of dimension  with


where the matrix components  contain the known coefficients . Equation (6) can now be rewritten as


where the array  contains the components of the yet unknown vectors  sorted (by choice) in the following way


It is required that (8) is always satisfied irrespective of , which is the case for the right null vectors of the matrix . The right null vectors are determined with a SVD, see e.g. [11], yielding the set of null vectors , where the index  denotes the number of the null vector.

The number of null vectors of  depends on the degrees  and  of the basis  and , and the coefficients  of the specific curve considered. For subsequent computations the number of null vectors must be more than the number of intersections of the curve with a line (or its algebraic degree). The non-square matrix  has  rows and  columns. Hence, its number of right null vectors must be equal or greater than .111 The number of right null vectors is larger then, for instance, a quadratic curve is described with a cubic polynomial (i.e. but ). More precisely, if is the largest integer such that then the number of right null vectors of  is exactly . In order to obtain the  intersections it is necessary to have


The set of null vector denoted with  introduced in (3) yields a set of moving lines


with a common intersection point on the curve . As an example, in Figure (a)a the description of a quadratic curve with  by two moving lines  and  with  is shown.

Next, the intersection of a given parametric line


with the curve  is considered, where and are two given vectors. We require that the determined set moving lines  and  intersect at the same point to be yet determined. According to (5) at the common intersection point it is required that


These equations describe the intersection of each moving line  with the given line  and can be rewritten as


As discussed the number of moving lines satisfies  and all of them can be combined in one homogenous equation system


with  and  representing the components of two matrices  and . Both matrices have  rows and  or more columns. To obtain  and  that satisfy (15) the following generalised eigenvalue problem is considered


That is, the eigenvalues  are the parametric coordinates of the intersection points on the line 

and the eigenvectors are (up to a multiplicative constant) the basis functions 

evaluated at the intersection points . Unfortunately, the matrices  and  are not always square and computing the values that satisfy this equation requires non-standard linear algebra techniques. However, for the purposes of intersection computation it is sufficient to consider a square eigenvalue problem obtained from (16) by taking only some of its columns. The non-complex eigenvalues of this square eigenvalue problem contain all the intersection points between the given line and the curve. As will be discussed further below, some of these non-complex eigenvalues may not be actual intersection points, but they can easily be identified.

A non-complex eigenvalue of (16), or its corresponding square eigenvalue problem, gives the potential intersection point . The respective unknown parameter value on the curve satisfies the equation , which is a hard to solve nonlinear problem. According to (15) and (16), however, the left null vector  corresponding  is proportional to the vector , which is exploited to determine . More precisely, if there is a single parameter value the ratio of any two consecutive components yields


It is assumed here that the monomials in  are sorted in increasing order so that the ratio of two consecutive entries is simply . If there are  parameter values corresponding to the eigenvalue , then the corresponding left null vector space is given by a matrix with rows. Keeping the previous assumption on the ordering of the monomials in , we define the matrices by taking the columns to of . Then, the parameter values are obtained by solving a generalised eigenvalue problem


2.2 Illustrative example

The intersection of a cubic Lagrange curve with a line depicted in Figure 6 is considered next. The cubic curve interpolating the points , , and is expressed in power basis with


The intersections between the given line


and the curve  are sought.

(a) Intersection points given by (26)
(b) Intersection points given by (27)
Figure 6: Intersection points between a cubic Lagrange curve and a line.

As discussed, the degree of the auxiliary polynomial needs to satisfy . In this example choosing  yields matrices and that have the dimensions and it is straightforward to compute the eigenvalues  of the generalised eigenvalue problem (16). To illustrate the more challenging case with non-square matrices  and , which turns out to be inevitable in the case of surfaces, we choose here  such that


The right null vectors  that define the moving lines are obtained from (8) with  and the components


This matrix has the dimensions , its rank is seven and the dimension of its right null space is five. Each of the corresponding five independent null vectors defines one moving line. Figure (b)b depicts the first four moving lines , , and at the parameter values and .

Substituting the line (20) as in (15) gives the two non-square matrices (with four significant digits)




One approach to obtaining the generalised eigenvalues of matrices and is to use pencil reduction, see [12, 13, 14]

, which is not a widely used linear algebra operation and may introduce additional numerical issues because of several numerical rank estimations. Alternatively, the eigenvalue problems defined by square submatrices

and of the largest size, e.g. the first four columns of  and , can be considered


Although there can only be three intersection points for a cubic curve this problem has four eigenvalues and eigenvectors. Notice that each column of , i.e , represents the intersection between the line  and the moving line . Taking a different couple of submatrices and of largest size, like the last four columns of  and , yields a different set of intersection points. Both sets of intersection points contain the three true intersection points in addition to one fictitious intersection point.

To demonstrate this observation, the intersection points with two different pairs of largest square submatrices and are computed. The eigenvalues for the problem defined by the first four columns of  and  are


and the eigenvalues for the problem defined by the last four columns  and  are


It is evident that for the first problem and for the second problem correspond to fictitious intersection points while the other three eigenvalues correspond to the true intersection points, see Figure 6. The coordinates of the three true intersection points are computed by introducing  in the line equation (20) yielding


The fictitious points can also be detected without computing several eigenvalue problems and comparing their eigenvalues. This is accomplished by determining the parametric coordinates of the intersection points on the curve . As indicated in (17) the parametric coordinates are computed using the eigenvectors. As an example, consider the eigenvalue problem (25) defined by first four colums of  and . Its eigenvalues are given on (26) and the coordinates of the intersection points are


The respective parameters of the intersection points are


A point is a true intersection point if and only if ; otherwise, it is a fictitious point. It can easily be found such that the fourth point is not an intersection point. Thus, in practice the intersection points are computed from a single couple of square matrices and .

3 Intersection of lines with surfaces

The extension of the introduced method to surfaces is straightforward. Let be a parametric surface, with , of bi-degree given either in Lagrange basis or power basis  with


where and are multi-indices, are the parametric surface coordinates, and and are the coefficients in the two basis. Usually, in finite element applications the degrees  and  of the surface  are the same.

In line with the curve case, a point on the surface is defined as the intersection of several moving planes of the form


where is an auxiliary vector collecting the parameters of the plane. Although  is now a plane instead of a line, it is still denoted with the same symbol to keep the notation simple. The parameters  are assumed to be of the following form


The bi-degree  of the power basis  has to be sufficiently high to describe all the intersection points of the surface  with a line.

The planes describing the surface  have to satisfy


which can be rearranged to


with a new power basis of bi-degree  and the array containing the sorted components of the vectors . The matrix has rows and columns. It has at least right null vectors, i.e. the difference between the number of columns and rows. The surface  has the algebraic degree , which is equal to its number of intersections with a line. In order to obtain all the intersections the following condition has to be satisfied


As in the curve case, the degree of  along the direction can be chosen as  yielding


By symmetry, it is also valid to choose and . After the right null vectors of  are computed the subsequent steps in computing the intersections are identical to the curve case.

Finally, for the sake of completeness, we mention that the introduced intersection algorithm also applies to triangular finite elements. If is a triangular parametric surface of degree then the degree of the auxiliary vector has to be chosen to satisfy ; see also [9, §3].

4 Conclusions

The introduced method is able to determine in one-shot all the intersection points between a line and a curve or surface. It can be applied to curves or surfaces given in any polynomial basis, like the Lagrange or Bernstein, after a straightforward conversion to the power basis. Aspects requiring further research include cases for which the eigenvalue problem (18) is degenerate and the preconditioning of the eigenvalue problem (25). In the accompanying implementation several choices have been considered for both.


  • [1] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 2005; 194:4135–4195.
  • [2] Karniadakis G, Sherwin S. Spectral/hp element methods for computational fluid dynamics. Oxford University Press, 2013.
  • [3] Wriggers P, Zavarise G. Computational contact mechanics. Encyclopedia of computational mechanics 2004; .
  • [4] Xie ZQ, Sevilla R, Hassan O, Morgan K. The generation of arbitrary order curved meshes for 3d finite element analysis. Computational Mechanics 2013; 51:361–374.
  • [5] Turner M, Peiró J, Moxey D. Curvilinear mesh generation using a variational framework. Computer-Aided Design 2018; 103:73–91.
  • [6] Rüberg T, Cirak F. Subdivision-stabilised immersed b-spline finite elements for moving boundary flows. Computer Methods in Applied Mechanics and Engineering 2012; 209–212:266–283.
  • [7] Patrikalakis NM, Maekawa T. Shape interrogation for computer aided design and manufacturing. Springer, 2009.
  • [8] Sederberg TW. Computer Aided Geometric Design. Computer Aided Geometric Design Course Notes, Brigham Young University, 2012.
  • [9] Busé L. Implicit matrix representations of rational Bézier curves and surfaces. Computer-Aided Design 2014; 46:14–24.
  • [10] Sederberg TW, Chen F. Implicitization using moving curves and surfaces. SIGGRAPH 1995 Conference Proceedings, ACM, 1995; 301–308.
  • [11] Strang G. Introduction to Linear Algebra. Fifth edn., Wellesley-Cambridge Press, 2016.
  • [12] Xiao X, Sabin M, Cirak F. Interrogation of spline surfaces with application to isogeometric design and analysis of lattice-skin structures. Submitted for publication. 2018; .
  • [13] Busé L, Luu Ba T. Matrix-based Implicit Representations of Rational Algebraic Curves and Applications. Computer Aided Geometric Design 2010; 27(9):681–699.
  • [14] Beelen T, Van Dooren P. An improved algorithm for the computation of Kronecker’s canonical form of a singular pencil. Linear Algebra and its Applications 1988; 105:9–65.