surface_intersection module¶

This is a collection of procedures and types for computing intersections between two Bézier surfaces in \(\mathbf{R}^2\). The region(s) of intersection (if non-empty) will be curved polygons defined by the Bézier curve segments that form the exterior.

Note

In most of the procedures below both the number of nodes \(N\) and the degree \(d\) of a Bézier surface are provided. It is redundant to require both as arguments since \(N = \binom{d + 2}{2}\). However, both are provided as arguments to avoid unnecessary re-computation, i.e. we expect the caller to know both \(N\) and \(d\).

Procedures¶

void locate_point_surface(int *num_nodes, double *nodes, int *degree, double *x_val, double *y_val, double *s_val, double *t_val)¶

This solves the inverse problem \(B(s, t) = (x, y)\) (if it can be solved). Does so by subdividing the surface until the sub-surfaces are sufficiently small, then using Newton’s method to narrow in on the pre-image of the point.

This assumes the surface is “valid”, i.e. it has positive Jacobian throughout the unit triangle. (If this were not true, then multiple inputs could map to the same output.)

Parameters

num_nodes (int*) – [Input] The number of nodes \(N\) in the control net of the Bézier surface.
nodes (double*) – [Input] The actual control net of the Bézier surface as a \(2 \times N\) array. This should be laid out in Fortran order, with \(2 N\) total values.
degree (int*) – [Input] The degree \(d\) of the Bézier surface.
x_val (double*) – [Input] The \(x\)-value of the point being located.
y_val (double*) – [Input] The \(y\)-value of the point being located.
s_val (double*) – [Output] The first parameter \(s\) of the solution. If \((x, y)\) can’t be located on the surface, then s_val = -1.0.
t_val (double*) – [Output] The second parameter \(t\) of the solution. If \((x, y)\) can’t be located on the surface, then this value is undefined.

Signature:

void
locate_point_surface(int *num_nodes,
                     double *nodes,
                     int *degree,
                     double *x_val,
                     double *y_val,
                     double *s_val,
                     double *t_val);

void newton_refine_surface(int *num_nodes, double *nodes, int *degree, double *x_val, double *y_val, double *s, double *t, double *updated_s, double *updated_t)¶

This refines a solution to \(B(s, t) = (x, y) = p\) using Newton’s method. Given a current approximation \((s_n, t_n)\) for a solution, this produces the updated approximation via

\[\begin{split}\left[\begin{array}{c} s_{n + 1} \\ t_{n + 1} \end{array}\right] = \left[\begin{array}{c} s_n \\ t_n \end{array}\right] - DB(s_n, t_n)^{-1} \left(B(s_n, t_n) - p\right).\end{split}\]

Parameters

num_nodes (int*) – [Input] The number of nodes \(N\) in the control net of the Bézier surface.
nodes (double*) – [Input] The actual control net of the Bézier surface as a \(2 \times N\) array. This should be laid out in Fortran order, with \(2 N\) total values.
degree (int*) – [Input] The degree \(d\) of the Bézier surface.
x_val (double*) – [Input] The \(x\)-value of the point \(p\).
y_val (double*) – [Input] The \(y\)-value of the point \(p\).
s (double*) – [Input] The first parameter \(s_n\) of the current approximation of a solution.
t (double*) – [Input] The second parameter \(t_n\) of the current approximation of a solution.
updated_s (double*) – [Output] The first parameter \(s_{n + 1}\) of the updated approximation.
updated_t (double*) – [Output] The second parameter \(t_{n + 1}\) of the updated approximation.

Signature:

void
newton_refine_surface(int *num_nodes,
                      double *nodes,
                      int *degree,
                      double *x_val,
                      double *y_val,
                      double *s,
                      double *t,
                      double *updated_s,
                      double *updated_t);

void surface_intersections(int *num_nodes1, double *nodes1, int *degree1, int *num_nodes2, double *nodes2, int *degree2, int *segment_ends_size, int *segment_ends, int *segments_size, CurvedPolygonSegment *segments, int *num_intersected, SurfaceContained *contained, Status *status)¶

Compute the intersection of two Bézier surfaces. This will first compute all intersection points between edges of the first and second surface (nine edge pairs in total). Then, it will classify each point according to which surface is “interior” at that point. Finally, it will form a loop of intersection points using the classifications until all intersections have been used or discarded.

Tip

If the status returned is INSUFFICIENT_SPACE that means either

segment_ends_size is smaller than num_intersected so segment_ends needs to be resized to at least as large as num_intersected.
segments_size is smaller than the number of segments. The number of segments will be the last index in the list of edge indices: segment_ends[num_intersected - 1]. In this case segments needs to be resized.

This means a successful invocation of this procedure may take three attempts. To avoid false starts occurring on a regular basis, keep a static workspace around that will continue to grow as resizing is needed, but will never shrink.

Parameters

num_nodes1 (int*) – [Input] The number of nodes \(N_1\) in the control net of the first Bézier surface.
nodes1 (double*) – [Input] The actual control net of the first Bézier surface as a \(2 \times N_1\) array. This should be laid out in Fortran order, with \(2 N_1\) total values.
degree1 (int*) – [Input] The degree \(d_1\) of the first Bézier surface.
num_nodes2 (int*) – [Input] The number of nodes \(N_2\) in the control net of the second Bézier surface.
nodes2 (double*) – [Input] The actual control net of the second Bézier surface as a \(2 \times N_2\) array. This should be laid out in Fortran order, with \(2 N_2\) total values.
degree2 (int*) – [Input] The degree \(d_2\) of the second Bézier surface.
segment_ends_size (int*) – [Input] The size of segment_ends, which must be pre-allocated by the caller.
segment_ends (int*) – [Output] An array (pre-allocated by the caller) of the end indices for each group of segments in segments. For example, if the surfaces intersect in two distinct curved polygons, the first of which has four sides and the second of which has three, then the first two values in segment_ends will be [4, 7] and num_intersected will be 2.
segments_size (int*) – [Input] The size of segments, which must be pre-allocated by the caller.
segments (CurvedPolygonSegment*) – [Output] An array (pre-allocated by the caller) of the edge segments that make up the boundary of the curved polygon(s) that form the intersection of the two surfaces.
num_intersected (int*) – [Output] The number of curved polygons in the intersection of two surfaces.
contained (SurfaceContained*) – [Output] Enum indicating if one surface is fully contained in the other.
status (Status*) –
[Output] The status code for the procedure. Will be
- SUCCESS on success.
- INSUFFICIENT_SPACE if segment_ends_size is smaller than num_intersected or if segments_size is smaller than the number of segments.
- UNKNOWN if the intersection points are classified in an unexpected way (e.g. if there is both an ignored corner and a tangent intersection, but no other types).
- NO_CONVERGE if the two curves in an edge pair don’t converge to approximately linear after being subdivided 20 times. (This error will occur via curve_intersections().)
- An integer \(N_C \geq 64\) to indicate that there were \(N_C\) pairs of candidate segments during edge-edge intersection that had overlapping convex hulls. This is a sign of either round-off error in detecting that the edges are coincident curve segments on the same algebraic curve or that the intersection is a non-simple root. (This error will occur via curve_intersections().)
- BAD_MULTIPLICITY if the two curves in an edge pair have an intersection that doesn’t converge to either a simple or double root via Newton’s method. (This error will occur via curve_intersections().)
- EDGE_END If there is an attempt to add an intersection point with either the \(s\) or \(t\)-parameter equal to 1 (i.e. if the intersection is at the end of an edge). This should not occur because such intersections are “rotated” to the beginning of the neighboring edge before the boundary of the curved polygon is formed.
- SAME_CURVATURE if the two curves in an edge pair have identical curvature at a tangent intersection.
- BAD_INTERIOR if a curved polygon requires more than 10 sides. This could be due to either a particular complex intersection, a programming error or round-off which causes an infinite loop of intersection points to be added without wrapping around back to the first intersection point.

Signature:

void
surface_intersections(int *num_nodes1,
                      double *nodes1,
                      int *degree1,
                      int *num_nodes2,
                      double *nodes2,
                      int *degree2,
                      int *segment_ends_size,
                      int *segment_ends,
                      int *segments_size,
                      CurvedPolygonSegment *segments,
                      int *num_intersected,
                      SurfaceContained *contained,
                      Status *status);

void free_surface_intersections_workspace(void)¶

This frees any long-lived workspace(s) used by libbezier throughout the life of a program. It should be called during clean-up for any code which invokes surface_intersections().

Signature:

void
free_surface_intersections_workspace(void);

Types¶

CurvedPolygonSegment¶

Describes an edge of a CurvedPolygon formed when intersecting two curved Bézier surfaces. The edges of the intersection need not be an entire edge of one of the surfaces. For example, an edge \(E(s)\) may be restricted to \(E\left(\left[\frac{1}{4}, \frac{7}{8}\right]\right)\).

double start¶: The start parameter of the segment. In the restriction \(E\left(\left[\frac{1}{4}, \frac{7}{8}\right]\right)\), the start would be 0.25.

double end¶: The end parameter of the segment. In the restriction \(E\left(\left[\frac{1}{4}, \frac{7}{8}\right]\right)\), the end would be 0.875.

int edge_index¶: An index describing which edge the segment falls on. The edges of the first surface in the intersection are given index values of 1, 2 and 3 while those of the second surface are 4, 5 and 6.

In the header bezier/surface_intersection.h, this is defined as

typedef struct CurvedPolygonSegment {
  double start;
  double end;
  int edge_index;
} CurvedPolygonSegment;

SurfaceContained¶

This enum is used to indicate if one surface is contained in another when doing surface-surface intersection.

NEITHER¶: (0) Indicates that neither surface is contained in the other. This could mean the surfaces are disjoint or that they intersect in a way other than full containment.

FIRST¶: (1) Indicates that the first surface (arguments will be ordered) is fully contained in the second. This allows for points of tangency, shared corners or shared segments along an edge.

SECOND¶: (2) Indicates that the second surface (arguments will be ordered) is fully contained in the first. This allows for points of tangency, shared corners or shared segments along an edge.