6
$\begingroup$

I’ve been racking my brain all afternoon on trying to solve this…

I have a set of noisy co-ordinate data

data = {{x1,y1,z1}…..{xn,yn,zn}}

that forms an ellipsoid . I’d like to transform the co-ordinate data into a unit sphere. As an illustration, my original data is in blue, and I’d like to transform it to the red data.

transform from blue to red: What is the best way to find the transformation coefficients

affine = {{sxx,sxy,sxz}, {syx, syy,syz}, {szx,szy,szz}} 
offset = Table[{ox,oy,oz}, {i,1,Length[data]}]  

such that:

transform = affine.data + offset 

subject to the constraint that

x^2+y^2+z^2 = 1

I assume this requires some kind of least square fit. I can’t figure out how I can get *fit * or NMinimize to make this work. Maybe I am using the wrong functions….

Ideas greatly appreciated….

Thanks, P

$\endgroup$
4
  • $\begingroup$ Physically, I know there are rotations involved. The data is essentially gravity along x,y,z axes… $\endgroup$
    – Pam
    Aug 15, 2013 at 21:41
  • $\begingroup$ So in general case you may want to adapt this method and with model of ellipsoid you have the transformation $\endgroup$
    – Kuba
    Aug 15, 2013 at 21:47
  • $\begingroup$ Don't forget to upvote any helpful answers...keeps everyone motivated... :) $\endgroup$
    – cormullion
    Aug 16, 2013 at 10:31
  • $\begingroup$ Also (saw this too late to answer) you may want to look at SingularValueDecomposition. $\endgroup$
    – gpap
    Aug 16, 2013 at 10:33

2 Answers 2

8
$\begingroup$

The obvious approach of forming a least-squares problem and throwing FindMinimum at it totally works; it just needs a little bit of a nudge in the right direction.

First, let me generate some points distributed on an ellipsoid, with some noise added to make it interesting:

sphere = Normalize /@ RandomVariate[NormalDistribution[], {1000, 3}];
a = RandomVariate[NormalDistribution[], {3, 3}];
b = RandomVariate[NormalDistribution[], 3];
data = a.# + b + 0.1 RandomVariate[NormalDistribution[], 3] & /@ sphere;
ListPointPlot3D[data, BoxRatios -> Automatic]

enter image description here

Here's your transformation, lightly modified:

affine = {{sxx, sxy, sxz}, {syx, syy, syz}, {szx, szy, szz}};
offset = {ox, oy, oz};
transform = affine.# + offset & /@ data;

Now let's just take the criterion $\|\mathbf x_i\|^2=1$, define its total squared error $\sum_i \big(\|\mathbf x_i\|^2-1\big)^2$, and see how well FindMinimum can minimize it:

sol = FindMinimum[Total[(Norm[#]^2 - 1)^2 & /@ transform], Flatten[{affine, offset}]];
ListPointPlot3D[transform /. sol[[2]], BoxRatios -> Automatic]

enter image description here

Oops. What happened? The minimization has shrunk all the points down to a tiny patch on the surface of the unit sphere. Well, that does make all the points have $\|\mathbf x_i\|^2\approx1$, but that's not really what we want. What this shows is that our objective has multiple local minima, and if we want to get the right one we should pick a better starting point. A good first guess for the offset is probably the one that makes the data set centered around zero:

min = Min /@ Transpose[data];
max = Max /@ Transpose[data];
mid = (min + max)/2;

Then we try again knowing a good starting point for the variables, and voilà:

sol2 = FindMinimum[Total[(Norm[#]^2 - 1)^2 & /@ transform], 
   Transpose[{Flatten[{affine, offset}], Flatten[{IdentityMatrix[3], -mid}]}]];
ListPointPlot3D[transform /. sol2[[2]], BoxRatios -> Automatic]

enter image description here

$\endgroup$
4
  • $\begingroup$ brilliant. The key was nudging it in direction of the offset… $\endgroup$
    – Pam
    Aug 16, 2013 at 10:34
  • $\begingroup$ Still getting used to mathematica, any way of figuring out what the best fit transform (i.e. affine and offset) values are. Would help me tune out my experimental procedure… $\endgroup$
    – Pam
    Aug 16, 2013 at 11:02
  • $\begingroup$ Take a look at the value of sol2; it is a pair containing {minimal value of objective function, rules for values of optimization variables at minimum}. So you can get the value of sxx out by evaluating sxx /. sol2[[2]], and so on. Or simply do affine /. sol2[[2]] and offset /. sol2[[2]]. $\endgroup$
    – user484
    Aug 16, 2013 at 11:08
  • $\begingroup$ tnx… that was quick.. ! $\endgroup$
    – Pam
    Aug 16, 2013 at 11:09
4
$\begingroup$

I'm not completely sure if this is the type of transform that you are looking for, but the following might help at least with fitting a set of data to an ellipsoid. First, generate some data. I'm using the parametric equations for an ellipsoid just because I just learned them ;-).

data = With[{a = 1, b = 2, c = 1.5}, 
    Table[{a Cos[u] Sin[v], b Sin[u] Sin[v], c Cos[v]}, {u, 0, 2 Pi, 
      2 Pi/20}, {v, -2 Pi, Pi, Pi/30}]]~Flatten~{1, 2};
ListPointPlot3D[data, SphericalRegion -> True, PlotStyle -> Blue]

Fitting the data to the cartesian coordinate equation for an ellipsoid is straightforward, but we first need to append the value of the function (which in this case is 1) to each set of points. I do this to the data within the nonlinearmodelfit function:

nlm = NonlinearModelFit[
   Partition[Flatten@Riffle[data, 1], 4],
   x^2/a^2 + y^2/b^2 + z^2/c^2,
   {a, b, c},
   {x, y, z}];
nlm["BestFitParameters"]
(* {a -> 1., b -> 2., c -> 1.5} *)

Here's were I don't know what you are looking for. I presume that if you want to transform the data, you want to normalize the three coefficients a, b, c such that they are equal. (Note the first line of code below is very slow.)

newdata = # /({a, b, c} /. nlm["BestFitParameters"]) & /@ data;
nlm2 = NonlinearModelFit[
   Partition[Flatten@Riffle[newdata, 1], 4],
   x^2/a^2 + y^2/b^2 + z^2/c^2,
   {a, b, c},
   {x, y, z}];
nlm2["BestFitParameters"]
(* {a -> 1., b -> 1., c -> 1.} *)
ListPointPlot3D[{data, newdata}, SphericalRegion -> True, 
 PlotStyle -> {Blue, Red}]

enter image description here

I haven't played with the AspectRatio in this image, so my sphere (red) looks squished.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.