% [d,w,rx,cv,wx] = best_dir( x, [a , sx ] )
%
% Some points x, are observed and one assumes that they belong to
% parallel planes. There is an unknown direction d s.t. for each
% point x(i,:), one has :
%
% x(i,:)*d == w(j(i)) + noise
%
% where j is known(given by the matrix a ), but w is unknown.
%
% Under the assumption that the error on x are i.i.d. gaussian,
% best_dir returns the maximum likelihood estimate of d and w.
%
% This function is slower when cv is returned.
%
% INPUT :
% -------
% x : D x P P points. Each one is the sum of a point that belongs
% to a plane and a noise term.
%
% a : P x W 0-1 matrix describing association of points (rows of
% x) to planes :
%
% a(p,i) == 1 iff point x(p,:) belongs to the i'th plane.
%
% Default is ones(P,1)
%
% sx : P x 1 Covariance of x(i,:) is sx(i)*eye(D).
% Default is ones(P,1)
% OUTPUT :
% --------
% d : D x 1 All the planes have the same normal, d. d has unit
% norm.
%
% w : W x 1 The i'th plane is { y | y*d = w(i) }.
%
% rx : P x 1 Residuals of projection of points to corresponding plane.
%
%
% Assuming that the covariance of x (i.e. sx) was known
% only up to a scale factor, an estimate of the
% covariance of x and [w;d] are
%
% sx * mean(rx.^2)/mean(sx) and
% cv * mean(rx.^2)/mean(sx), respectively.
%
% cv : (D+W)x(D+W)
% Covariance of the estimator at [d,w] ( assuming that
% diag(covariance(vec(x))) == sx ).
%
% wx : (D+W)x(D*P)
% Derivatives of [w;d] wrt to x.
%
% Author : Etienne Grossmann <etienne@isr.ist.utl.pt>
% Created : March 2000
%