Re: Intersection line and its error of two best fitting planes

[lynx . abraxas]

On 15/05/10 16:12:04, Leo Butler wrote:
>
>
> On Sat, 15 May 2010, address@hidden wrote:
>
> < Dear Leo!
> <
> <
> < Many thanks for Your reply.
> <
>
> < Also I wonder what is a? It's not defined in Your code.

I did

a=n;

to get Your code to run.

Now I'm trying to understand it...

- The value of c is never used, is it?
c=p'*a;

- But c is the shortest distance of the plane to the origin?

-  vector  x  would be loaded with my measured x,y,z-point coordinates of each
plane?

- for ols: 1/norm(w) is the implied estimation of c?

- for sqp it is the fourth entry in the solution vector?

- why -1 in t=norm(alpha(1:3))-1 ???


> < Is the resulting plane the same  for  both  methods  and  only  the  error
differ?
>  < Looking  at  the nice visualization of Stefans solution, I'd say it looks
well
> < fitted!
>
> No, least squares estimates the normal vector, but the constant is only
> implied (it's the reciprocal of the length of the normal vector). The
> bulk of the error occurs in the estimate of the constant (not
> surprising, since the OLS objective function does not utilise the
> constant.)
>
> Thus, you are getting a plane close to parallel to the real plane, but
> you can do better by reducing the error in the constant.

So sqp kind of does an additional optimization of the ols result by  adjusting
the displacement constant?


>  <  Still  I'm not sure how the error of the planes (what ever method  used)
would
> < go  into   the   error   of   the  orientation  and  displacement  of  the
intersection
> < line.
>
> Nonlinearly. As I said, this is a (differential) geometric problem.
>
> I would go about it somewhat differently.

OK. I'm trying to follow that although it's much less intuitive for me:-(

>
>  Let's think of the following map f: given two non-parallel affine planes in
R^3,
> P & Q, the map gives you P intersected with Q = an affine line L in R^3.
>
> You are given noisy data P' and Q' and want to estimate f(P,Q)=L by L'.

> The most natural thing to do is to set this up as a constrained least
> squares problem. First, an affine line L in R^3 has a unique point p that
> is closest to 0, and a unit direction vector v that is unique up +-1,
> so L={p+t*v : t in R}. Second, as you note, 3 points uniquely determine
> each affine plane and (assuming they are not parallel), an affine line.
> Third, these affine lines will not agree, but they can't be averaged
> (it's non-sense), so let's write an objective function that makes sense:

I can follow You till here but the function G down there I don't understand.

> G(L;{L_i}) = sum_{i,j} |p_i-p|^2 + d(v_i,v)^2
>
> where d is the smallest of v_i-v and v_i+v and L_i (L) is the i-th
> determined by p_i and v_i (p and v).


Well I think it all comes down on how to define my obj-function in octave.  So
with  x1  containing  my measured xyz-coordinates of points in plane P1 and x2
the points of the other plane P2 I want a function that optimzes my  direction
vector v and displacement vector p of the intersection line L. So I guess I've
to split up alpha into v and p.
But where do to get n1,c1 and n2,c2 from? From two previous runs  of  sqp  for
just the plane parameters???
And then I'm not sure how to calculate t in this case...


function t=obj(alpha)
  global x1 x2 N;
  p=alpha(4:6);
  v=alpha(1:3);
  v=v/norm(v);
  t=(x1*n1-c1)' * (x1*n1-c1)/N ...???
endfunction



>
> With your map f, you get an objective function that is a function of
> your data, that can be minimised to get L (i.e. p and v).


So if I had my obj function all right I'd just run
v=sqp([1;0;0;1;0;0],@obj,@norm_constraint,[])

And the result would be [v,p] of my line L?


> Finally, can I ask what you are using this for?

I'm  trying  to  determin  the  rotation  axis  of  a goniometer of a TEM. The
goniometer can translate the sample in x,y,z direction and tilt it  around  an
axis  close  to  y.  But  only  close.  This is why I'm trying to determin the
rotation axis with points from a flat sample that I measured at 45°  and  -45°
tilt.  Actually  this problem was ment to be solved within 1 or 2 hours but it
seems much more complicated than I thought:-(


Many thanks for Your help
Lynx