Singular Curves on Surfaces Algorithm (Sing3D)

OVERVIEW:

This algorithm traces singular curves on the level set of an arbitrary smooth function f: R3 -> R. It also tries to fill in the local structure of the level set in a neighborhood of the singular curves.

BASIC EXPLANATION OF ALGORITHM

A. COMPUTING THE SINGULAR CURVE

Figure 1: The surface is defined implicitly by (x2 + y2+z2-0.5) *(z-0.25)=0, which is the union of a sphere and a plane. It has a singular curve, the circle. The purple segments are tangents to the singular curve at the initial points from which we trace the curve.

The basic idea is essentially the same as the standard predictor corrector algorithm:
1) find one point on the curve;
2) compute a tangent vector to the curve at that point;
3) step out a small amount along the tangent vector;
4) relocate the curve at a new point, and store the singular data for later use;
5) go to step (2).

We want to emphasize that the singular curve (if it exists) is the solution of the over-determined system of the equations:

f=0,
fx=0,
fy=0,
fz=0.
Let F be the function F=(fx,fy,fz), then DF=hessian of f. If x0 is a solution of the system above and if DF(x0) has full rank, then x0 is an isolated solution. So if there is a solution curve passing through x0, DF(x0) must have rank less than 3. In the case that rank(DF(x0))=2, the tangent space to the singular curve at x0 coincides with the null space of DF(x0). (We actually use this fact to process step (2).) Currently we assume that rank(DF(x0))=2 except at finitely many points.

In order to accomplish step (1), we first use Newton's method on the system formed by F to find roots of F. Each root is a candidate solution to (f,F)=(0,0). As we attempt to iterate onto the singular curve, Newton's method will eventually complain at some point, say xk, since DF has rank less than or equal 2 along the singular curve. In this case, we restrict our searching for a candidate on the plane passing through the point xk and orthogonal to the null space of DF(xk). After we find a candidate, we then check to see if we really found a point on the singular curve or not. That is, we verify that f=0 and F=0.

For instance, let f=x2 - y2 then F=(2x,-2y,0) and DF is the diagonal matrix with 2, -2, and 0 on the main diagonal. Note that DF has always rank 2. Newton's method can not apply to the system F=0 to find the solutions. However the first two rows of DF are linearly independent and (0,0,1) is a vector in the null space of DF. So with an initial guess (a,b,c) we try to solve the system

2x=0,
-2y=0,
(x-a,y-b,z-c).(0,0,1)=0, that is z-c =0.
Therefore we find the candidate (0,0,c). Finally, we verify that f=0 and F=0, which are true. So we found a point on the singular curve.

If we encounter a point x0 on the curve with rank(DF) less than or equal to 1, then we will compute the intersection curves of the zero surface f=0 and a sphere of suitable size centered at x0. We will do the same thing for isolated singular points. This will help us to understand the local structure near those singular points.

In step (2), as we mentioned above, we use the null space of DF(x0) to compute the tangent of the singular curve at x0.

For steps (3) and (4), we pick two equations from the system F=0, based on which two row vectors of DF(x0) are linearly independent. Then we add the third equation: (x - w0) . T0 = 0 to make up a non-degenerate system, where w0 = x0 + step_length * T0 and T0 is a unit tangent vector at x0. Now we use Newton's method again to relocate a new point onto the curve using the new system of equations (of course, we need to check to see if we really got a solution for the original system or not).

We remark that we save the singular curve as a list of arcs for later use. Each arc contains a list of points on the singular curve and tangent vectors at each point.

B. FILL IN THE LOCAL STRUCTURE ALONG THE SINGULAR CURVE

Figure 2: The blue circles are the boundary curves. They are actually the intersection of the surface with a torus around the singular curve.

To fill in the local structure, we do the following
1) form a cylinder around the singular curve;
2) trace the boundary curve (the intersection of the level set and the cylinder), and store the data for the regular two-parameter continuation to trace the rest of the level surface;
3) triangulate the surface between the boundary curve and the singular curve.

Assume that s(t) is a parameterization of the singular curve. Then the cylinder of radius r in step (1) may be characterized by

(x-s(t)).(x-s(t))-r2=0;
(x-s(t)).s'(t)=0;
where x is a point on the cylinder. So the boundary curve can be defined by
(x-s(t)).(x-s(t))-r2=0;
(x-s(t)).s'(t)=0;
f(x)=0.
Note that the domain of the system is four dimensional (including t) and the range is three dimensional. Recall also that we already stored the singular data, so we can calculate s(t) and s'(t) by linear or quadratic interpolation for any given t. This is the perfect situation to apply the standard predictor-corrector algorithm.

We trace the boundary curve "arc by arc" as we traverse the singular curve. For example, let x0 and x1 be two end points of an arc on the singular curve. Let T0 and T1 be unit vectors tangent to the singular curve at x0 and x1 respectively. Finally let S0 and S1 be the circles of radius r, centered at x0 and x1 and perpendicular to T0 and T1 respectively. Then we use the points that are on S0 and S1 and simultaneously are on the zero surface of f as the initial points to trace parts of the boundary curve. The computation of these points reduces to a one dimensional root finding problem. So we can use a "bisection method" to find roots. The rest of step (2) is standard.

CONTROLLING THE ALGORITHM FROM PISCES:

The Pisces panel to the sing3D algorithm currently has 11 parameters which may be set from the algorithm window:

The meaning of the five parameters below is pretty much the same as the standard predcorr algorithm. They were documented in detail in the documentation for the Predictor-Corrector Algorithm.

1. Step_Length Default value: 0.02
2. Max_Steps Default value: 120
3. Random_Points Default value: 0
4. Selected_Points Default value: 0
5. Uniform_Points Default value: 3

The following parameters are specific to Sing3D algorithm.

6. Display_IP Default value: 1
This field should contain a nonnegative integer value. If the value is positive, the algorithm display those initial points from which we start to trace the curve.

7. Radius Default value: 0.3
This field should contain a positive real value. It is the distance between the singular curve and the cylinder surrounding it.

8. Zero_Vector_Tol Default value: 0.001
This field should contain a positive real value. It is the tolerant error for identifying a zero vector. That is, if the absolute value of each component of a vector v is less than this field, then the vector is treated as the zero vector.

9. Boundary_Smoother Default value: 1.5
This field should contain a positive real value. It is the expanding factor to insure that the maximum number of steps taken to trace the boundary curve is appropriate so that the boundary curve is smoother. In general, the larger the value, the smoother the curve.

10. Circle_Division Default value: 100
This field should contain a nonnegative integer value. It is the number of divisions for the bisection method to find initial points on the circle orthogonal to tangent vectors at the end of arcs of the singular curve.

11. Trace_Boundary_Curve Default value: 1
This field should contain a nonnegative integer value. If the value is positive, the algorithm will plot the boundary curve otherwise it won't.

KNOWN BUGS:

BUG REPORTS:

software@geom.umn.edu

ACKNOWLEDGEMENTS:

This algorithm was jointly developed and implemented by Frederick J. Wicklin and Chia-Hsing Nien.
Chia-Hsing Nien<nien@geom.umn.edu>
Last modified: Wed Apr 24 19:19:54 1996