help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Leo Butler
Subject: Re: Intersection line and its error of two best fitting planes
Date: Fri, 14 May 2010 21:25:46 +0100 (BST)
User-agent: Alpine 1.10 (DEB 962 2008-03-14)


On Fri, 14 May 2010, Stefan Neumann wrote:

< 
< 
< 2010/5/14 <address@hidden>
<       Hello!
< 
< 
<       I have two sets of measured points (in 3D). Points of each set are 
supposed to
<       sit on a plane.
<       Now I want to find the intersection line and its  directional  and  
positional
<       error.
< 
< 
< The best fitting plane could be found with linear regression. The plane is in 
vector 'B', errors etc can be found in variables 'Residuals' and 'STATS'
 
It looks like a linear regression problem, but it's not.
An affine plane P in R^3 is determined by a unit normal n and
a constant c such that

x \in P iff <n,x> = c.

It's clear that (-n,-c) and (n,c) determine the same plane.
This says the set of planes in R^3 is S^2 \times R modulo +-1.

You can determine (n,c) up to sign by minimizing

f(n,c) = sum_i (<n,x_i> - c)^2

s.t.

<n,n>=1.

Code:

function t=obj(alpha)
  global x N;
  c=alpha(4);
  n=alpha(1:3);
  n=n/norm(n);
  t=(x*n-c)' * (x*n-c)/N;
endfunction

# alpha=[n;c]
function t=norm_constraint(alpha)
  t=norm(alpha(1:3))-1;
endfunction

## an example

global x N;
epsilon=1e-2;
N=10;
n=[1;2;-4];
v1=[2;-1;0];
v2=[0;2;1];
p=[0;0;0];
c=p'*a;
y=ones(N,1);
for i=1:N
  x(i,1:3)=p+rand(1)*v1+rand(1)*v2+epsilon*rand(1)*n;
endfor;
w=ols(y,x)                                #least squares
v=sqp([1;2;0;3],@obj,@norm_constraint,[]) #non-linear constrained
least-squares
norm(v(1:3)-a/norm(a))
norm(w/norm(w)-a/norm(a))

You should see that the error for the ols estimate is almost an
order of magnitude larger.

Leo
-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



reply via email to

[Prev in Thread] Current Thread [Next in Thread]