Re: fit plane to xz axis of data
- To: mathgroup at smc.vnet.net
- Subject: [mg89610] Re: fit plane to xz axis of data
- From: will parr <willpowers69 at hotmail.com>
- Date: Sat, 14 Jun 2008 05:32:12 -0400 (EDT)
> On Jun 11, 11:57 pm, "Szabolcs Horv=E1t"
> <szhor... at gmail.com> wrote:
> > On Wed, Jun 11, 2008 at 10:20, will parr
> <willpower... at hotmail.com> wrote:=
>
> > > Dear Math Forum,
> >
> > > I am having problems fitting a plane to some
> data. I am using the follow=
> ing procedure to fit the plane to my data (data is
> pasted at the bottom of t=
> his message):
> >
> > > In[2]:= plane = Fit[data, {1, x, y}, {x, y}]
> >
> > > Out[2]= 3.28723- 0.189001 x - 0.0874557 y
> >
> > > Then the following to display the points and
> plane:
> >
> > > Show[ListPointPlot3D[data, Boxed -> False, Axes
> -> True,
> > > AxesEdge -> {{-1, -1}, {-1, -1}, {-1, -1}},
> > > AxesLabel -> {"x", "y", "z"}],
> > > Plot3D[plane, {x, Min[data[[All, 1]]],
> Max[data[[All, 1]]]}, {y,
> > > Min[data[[All, 2]]], Max[data[[All, 2]]]},
> > > PlotStyle -> Directive[Green, Opacity[0.5]],
> > > Mesh -> None]]
> >
> > > This works fine at fitting the plane in the xy
> axis, but I want to fit t=
> he plane to the xz axis. Does anyone know how this is
> done?
> >
> > You just need to permute the coordinates:
> >
> > Replace[data, {x_, y_, z_} :> {x, z, y}, 1]
> >
> > Or, if necessary, you can do an "isotropic" fitting
> (that prefers no
> > axis) by calculating the distances from the plane.
> The original
> > fitting attempt showed that the plane will not pass
> through {0,0,0} so
> > we can use the equation a x + b y + c z == 1. The
> distance-squared of=
>
> > point {x,y,z} from this plane is (a x + b y + c z -
> 1)^2/(a^2 + b^2 +
> > c^2), so the fitting can be done like this:
> >
> > In[2]:= planei =
> > a x + b y + c z /.
> > Last@NMinimize[
> > Expand@Total[
> > Function[{x, y, z}, (a x + b y + c z - 1)^2]
> @@@ data]/(
> > a^2 + b^2 + c^2), {a, b, c}]
> >
> > Out[2]= 0.100048 x + 0.0136288 y + 0.0112116 z
> >
> > In[3]:= plane = Fit[data, {1, x, y}, {x, y}]
> >
> > Out[3]= 3.28723- 0.189001 x - 0.0874557 y
> >
> > In[4]:= plane2 =
> > Fit[Replace[data, {x_, y_, z_} -> {x, z, y}, 1],
> {1, x, z}, {x, z}]
> >
> > Out[4]= 8.12817- 0.817828 x - 0.225109 z
> >
> > In[5]:= {{xa, xb}, {ya, yb}, {za, zb}} = {Min[#],
> Max[#]} & /@
> > Transpose[data]
> >
> > Out[5]= {{7.08299, 11.4614}, {-5.66284, 4.59986},
> {-1.83404,
> > 4.49073}}
> >
> > In[6]:= Show[
> > ContourPlot3D[planei == 1, {x, xa, xb}, {y, ya,
> yb}, {z, za, zb},
> > ContourStyle -> Directive[Blue, Opacity[.5]],
> Mesh -> None],
> > Plot3D[plane, {x, xa, xb}, {y, ya, yb},
> > PlotStyle -> Directive[Green, Opacity[0.5]], Mesh
> -> None],
> > ContourPlot3D[plane2 == y, {x, xa, xb}, {y, ya,
> yb}, {z, za, zb},
> > ContourStyle -> Directive[Red, Opacity[.5]], Mesh
> -> None],
> > ListPointPlot3D[data, PlotStyle -> Black],
> > PlotRange -> All, AxesLabel -> {"x", "y", "z"}
> > ]
> >
> > Your points don't seem to lie on a plane, so the
> three approaches give
> > give different results.
>
> The following will find the best-fitting (i.e.,
> least-squares) plane:
>
> m = Mean[data]; {u,w,v] =
> = SingularValueDecomposition[#-m&/@data];
>
> v is a rotation matrix. Tr[w,List] (i.e., the
> diagonal elements of w)
> gives the norms of the projections of the points onto
> the new axes.
> They are in decreasing order, so
>
> ({x,y,z}-m).v[[All,3]] == 0
>
> gives the equation which the best-fitting plane
> satisfies.
>
I'm really interested in this solution, but when I to use this method, I get an answer that I don't really understand, and cannot display (in the form of a plane).
using the original data in this message:
In[122]:= m = Mean[data]; {u, w, v} =
SingularValueDecomposition[# - m & /@ data];
In[123]:= Tr[w, List]
Out[123]= {16.3822, 9.63356, 5.71769}
In[124]:= ({x, y, z} - m).v[[All, 3]] == 0
Out[124]=
0.984797 (-9.86988 + x) + 0.134151 (0.269052+ y) +
0.110358 (-1.44535 + z) == 0
I don't understand this output, it is not of the form:
ax + by + cz
Additionally, is it possible to compute a vector running parallel to the long side (long axis of the SVD plane) of the computed plane, originating at the mean (centroid) of the data, with a length of 1?
I have the following bits of code with respect to this:
(*normalizes a vector to unit length*)
normalizevec[v_] := v/Sqrt[v.v + 0.0];
(*takes equation of a plane (e.g.1+2x+3y) and returns normal vector*)
planeeq2normalvec[e_] := normalizevec[{D[e, x], D[e, y], -1}];
I then use VectorAngle in the following way to find the 2D angle (in degrees) between the calculated 3D vector (that i don't have in this instance) and the x axis:
XVector = {{1, 0, 0}, {1, 0, 0}};
AngleX =
VectorAngle[
XVector[[1]][[{1, 2}]], (CALCULATED3DVECTOR[[{1, 2}]])]*180/\[Pi];
Thanks for your help,
best wishes,
Will