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

Wastingtime,

I apologize for not sending enough information--I forgot to include the
points. Only three are necessary however since the function doesn't look
beyond the first-last range. But you seem to have grasped what was going
on.

Thanks very much for going into that with such detail! Unfortunately I
flunked Linear algebra and that was almost 20 years ago. Higher math just
isn't my thing, unfortunately. Is there anything that a more logic-oriented
person could do do make this code reliable?

One thing I notice is that there are only three points in the range to which
it is trying to fit the curve--a beginning, middle, and end point. I'm not
sure but that may have been the case every time it has freaked out on me. I
guess that the tangent vectors are pointing away from each other instead of
relatively toward each other? That is resulting in the behavior you
describe, of the curve shooting away, turning around, shooting away in the
other direction, then coming back.

I didn't quite understand your suggestion at the end? Point on the line
closest to the origin? I guess what I need is just a catch to tell if it is
shooting away like that, so if the alpha values are large compared to the
distance between the endpoints, say larger than 4 or 5 times the distance,
do you think that would work?

Also, I wonder if the tangents are pointing away from each other due to some
bad point data. Maybe I can filter the points more before applying the
curve fitting. And, just to be super hacky about it, I could just draw
straight lines between the points if it's gotten down to a three-point
range. This is just a simple graphics application, it doesn't have to be
exactly perfect.

Thanks again,

Bint




On 8/29/05 11:01 PM, in article
zSQQe.4506$FW1.1103@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx, "wasting time"
<wastingtime@xxxxxxxxxxxxxxx> wrote:

> "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