Re: Detail: problem with FitCurves.c code from Graphics Gems byPhilip J. Schneider

"Bint" <juliosity@xxxxxxx> wrote in message
news:BF39138A.11DF6%juliosity@xxxxxxxxxx

> The function below, called with these arguments, does not work. It
> calculates too-large values for the control-point vector offsets. But
> why,
> and what can be done about it?
>
> Point2 *d (x = 307.3595195323904, y = 337.7928716903604)
> Int first = 18, last = 20
> Double *uprime = 0.000, 0.5389914414568432, 1.000
>
> tHat1 (h = -0.3326114973176150, v = 0.9430639383690451)
> tHat2 (h = 0.3273949711899671, -0.9448875768256881)

The "d" array had better contain at least 21 points since the code
accesses d[first] and d[last]. What you posted is not enough for
a reader to test the code. Regardless, though, read on.

I have doubts about the robustness of the source code. The block
of code
if (det_C0_C1 == 0.0) {
det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
}
is suspect because of the exact comparison of a floating-point value
to zero.

The fitting method in the Graphics Gems article leads to a linear
system of two equations in the two unknowns alpha_1 and alpha_2.
The 2-by-2 matrix of coefficients A = [a[i][j]] depends only on the
u[i] parameter values and the dot product of the unit tangents T1
and T2. The coefficients are
a[0][0] = 3 * sum_{i=1}^n (1-u[i])^4 * u[i]^2
a[0][1] = 3 * Dot(T1,T2) * sum_{i=1}^n (1-u[i])^3 * u[i]^3
a[1][0] = a[0][1]
a[1][1] = 3 * sum_{i=1}^n (1-u[i])^2 * u[i]^4
If the matrix is nearly singular (determinant is nearly zero), then you
can run into numerical problems as you have seen. For your data,
T2 is nearly opposite in direction to T1. The matrix coefficients are
the following (computed with my own code)
a[0][0] = 0.039365968924115882
a[0][1] = -0.046024300345917729
a[1][1] = 0.053810460354085050
and the determinant of A is 6.4687760039164321e-008, so the
matrix is nearly singular.

The hack for handling the singularity (checking if det_C0_C1 is zero)
is not generally a good idea. I think a better approach in the
implementation is to trap the case when the matrix is nearly singular
and reduce the dimensionality of the problem.

The S(alpha_1,alpha_2) function mentioned in the article has the graph
of a (nonnegative) "paraboloid" when the matrix A is not singular. This
paraboloid has a unique global minimum, which is what the code you
posted is attempting to compute. When A is singular, the graph of S is
a (nonnegative) parabolic cylinder, so you get an entire line of pairs
(alpha_1,alpha_2) which minimize S. The question is which pair to
choose. All provide the same minimum S, so the Graphics Gems
article can no longer help you.

Let the pair (alpha_1,alpha_2) be written as the 2-by-1 vector X.
The linear system provided by the article is of the form A*X = B,
where B is a 2-by-1 vector. The matrix A is symmetric, so it can
be factored using an eigendecomposition to A = R*D*R^T, where
R is a rotation matrix and D is a diagonal matrix. The notation R^T
means take the transpose of R. The columns of R are eigenvectors
for A and the diagonal entries of D are the corresponding eigenvalues.
You can then write R*D*R^T*X = B, or equivalently D*Y = C,
where Y = R^T*X and C = R^T*B.

In the event A is singular (in this application), one of the eigenvalues is
zero. In your example, my test code produced
D = diagonal(6.9425543059581374e-007, 0.093175735022770351)
In this case you may as well assume d[0][0] = 0. Also, R = [r[i][j]] with
r[0][0] = -0.75994254054604304
r[0][1] = 0.64999025767193286
r[1][0] = -r[0][1]
r[1][1] = r[0][0]
The geometry of the situation implies that D*Y = C has infinitely many
solutions (the line of (alpha_1,alpha_2) pairs) and D is singular, so it
must be the case that C = (c[0],c[1]) with c[0] = 0. That is, the
first row of [D|C] has all zeros. The second row produces one equation,
d[1][1]*y[1] = c[1], in which case
(eqn. 1) r[1][0]*alpha_1 + r[1][1]*alpha_2 = y[1] = c[1]/d[1][1]
In fact, this is the line of points which minimize S.

Now to reduce the dimensionality. The choices for alpha_1 and
alpha_2 tell you where to place V1 and V2 on the tangent lines at
V0 and V3. The path of a point on the cubic Bezier curve is a straight
line segment (this is essentially the geometric interpretation for the
singularity of A). If you imagine choosing alpha_1 and alpha_2 very
large, a particle traveling along this segment will "zig zag", moving
in one direction a large distance, and then backtracking, with yet
another change in direction to reach the final control point. So you
probably want to restrict the sizes of alpha_1 and alpha_2. For
example, choose (alpha_1,alpha_2) to be the point on the line
defined by (eqn. 1) which is closest to the origin.


.

Relevant Pages
Loading