[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 
Created : March 2000