       Re: 3D data set

• To: mathgroup at smc.vnet.net
• Subject: [mg50999] Re: [mg50948] 3D data set
• From: "Melko, Mike" <Mike.Melko at northern.edu>
• Date: Fri, 1 Oct 2004 04:48:21 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

```-----Original Message-----
From: Giovanni Bellesia [mailto:giovanni.bellesia at ucd.ie]
To: mathgroup at smc.vnet.net
Subject: [mg50999] [mg50948] 3D data set

Dear all,

I have a general question regarding a topic which is not completely new
to the forum.
I have a 3D data set (from a Monte Carlo simulation) which are supposed
to lay approximately on a helix.
Does anybody knows a clear and efficient procedure to fit these points
to a regular, circular helix.

Thanks

Giovanni

---- Response ------

Giovanni,

Off the top of my head, I don't think your problem has a simple
solution, even though it's easy to state. What follows is a description
of what I think you have to do. I haven't written things in Mathematica
syntax because I'm really just outlining an algorithm.

1. The space of helices in R^3 is six dimensional. We can represent it
as follows. Let (x1,x2,x3) be coordinates in R^3, and let (e1, e2, e3)
denote the standard basis. Now let u1 be a unit vector representing the
direction of the axis of the helix. Complete u1 to an orthonormal basis
(u1, u2, u3) in such a way that (e1, e2, e3) -> (u1, u2, u3) is
orientation preserving. The a helix c(t) can then be parameterized as
follows:

c(t) = x u2 + y u3 + r Cos(w t) u2 + r Sin(w t) u3 + t u1

x, y describe the position of the axis in the plane orthogonal to u1, w
describes the speed at which the helix coils, r is the radius of the
cylinder containing it, and theta, phi are the Euler angles of u1.

2. Let p be a point in the data set. Then the half of the distance
squared of p to c(t) is given by E(t,p) := (1/2) <c(t)-p, c(t)-p>.
Differentiate to find the extrema. This gives you the equation <c'(t),
c(t)-p> = 0. Unfortunately, I don't think this can be solved explicitly
for t, because it involves sums of sines and cosines. Note also, that
there are infinitely many local minima. (Draw a picture, and consider
the intersection of the plane containing p and the axis of c(t) and
intersect it with c(t). I claim the intersection points are the local
minima.)

3. Having simplified E(t,p) as much as possible, the function you want
to minimize is

E(t,D) := Sum[E(t,p), p in D)

where D is your data set. E(t,D) potentially has lots of local extrema,
and I think you'll have to resort to something like gradient descent or
Levenberg-Marquardt to find an approximate minimum. If this sounds like
a reasonable approach to you, you may want to consult "Numerical Recipes
in C" for useful discussions on gradient descent, Levenberg-Marquardt,
and their use in fitting nonlinear models to data.

I forgot to describe how the dependance of the minimizing function
depends on the variables (x,y,r,w,theta,phi). Suppose you could
explicitly find the time parameter t0 so that E(t0,p) is an absolute
minimum. Since c depends on x, y, r, w, theta, phi, this gives us a
function F(p; x,y,r,w,theta,phi) for the distance of p to an arbitrary
helix, in other words

F(p; x,y,r,w,theta,phi) := E(p; c(t0; x,y,r,w,theta,phi))

where c(t; x,y,r,w,theta,phi) is the helix with all parameters listed.

Now define

F(x,y,r,w,theta,phi) := Sum[F(p; x,y,r,w,theta,phi), p in D]

The latter is what you want to minimize over the variables
(x,y,r,w,theta,phi) using Levenberg Marquardt, say. Not being able to
solve for t0 in <c'(t), c(t)-p> = 0 explicitly is a major stumbling
block because F(p; x,y,r,w,theta,phi) = E(t0,p) is then not explicitly
defined. I don't see how to resolve this yet.

Mike Melko

```

• Prev by Date: Re: Re: Re: Hyperbolic function identity
• Next by Date: Re: Re: Re: Hyperbolic function identity
• Previous by thread: Re: 3D data set
• Next by thread: Re: 3D data set