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. I read something about this in a message dated may 2004 by D.L. but I wasn't able to download the related files 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