[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.