lans_orient.m

模式识别工具包

%	lans_orient	- Orient covariances of a principal surface
%
%	[psurf] = lans_orient(psurf,options)
%
%	_____OUTPUT_____________________________________________________________
%	psurf	principal surface with re-orientated covariances(structure)
%		see lans_psurf.m
%
%	_____INPUT______________________________________________________________
%	psurf	current principal surface			(structure)
%		see lans_psurf.m
%	options	parameter string				(string)
%		attenuate	amount of clipping along gradient direction
%				0 < attenuate < D/Q
%				=1 implies no orientation (spherical)
%	
%	_____SEE ALSO___________________________________________________________
%	lans_psurfinit lans_psurfgrad
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters
%	- only checks '-attenuate' in options
%	- assumed that psurf.covartype is spherical
%	- .covartype unimportant
%		creates internal fields
%		.beta	inverse of each of the M isotropic variances
%		.Cinv
%		.Cidet2
%	- for symmetric matrix A with eigenpair [V,D]
%		A	= V*D*V'
%		inv(A)	= V*(1/D)*V'
%		det(A)	= prod(diag(D));
%
%
%	(C) 1999.10.26 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

%	_____TO DO______________________________________________________________
%	- replace lans_gramschmidt with qr! (faster)
%	- handle repeated knots


function [psurf] = lans_orient(psurf,options)

if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin<2
	options	= [];
end
%-----	parameters
algo		= paraget('-algo',options);
alpha		= paraget('-attenuate',options);
covarmax	= paraget('-covarmax',options);
covarmin	= paraget('-covarmin',options);
covarclip	= paraget('-covarclip',options);

M		= psurf.M;
D		= psurf.D;
Q		= psurf.Q;
%_____	Catch Error
if alpha<=0
	error('Attenuation factor must be >0');
elseif alpha>=D/Q
	error('Attenuation factor must be <D/Q');
end
%__________	current gradient w.r.t. manifold
dfdx		= lans_psurfgrad(psurf);			% D x Q x M 

%_______________________________________________________________________________
if algo>2		% For GTM and TIB
	
beta		= psurf.beta;
if length(beta)==1
	beta		= ones(1,psurf.M)*beta;
	psurf.beta	= beta;
end
%_____	initialize
psurf.covars	= zeros(D,D,M);

%	Uses the UNIT gradient vector as base
if alpha~=1
	psurf.Cinv	= zeros(D,D,M);
	psurf.Cidet2	= zeros(M,1);
	for m=1:M
		% constants c1, c2
		c1	= alpha/beta(m);
		c2	= (D-alpha*Q)/(beta(m)*(D-Q));

%		psurf.Cinv(:,:,m)	= eye(D)*beta(m)/attenuate;

		% Set of Q parallel col vectors in R^D;
		plgrad	= dfdx(:,:,m);

		% Set of D parallel unit col eigenvectors
		eigvec	= lans_gramschmidt(plgrad);
		eigval	= diag([c1*ones(1,Q) c2*ones(1,D-Q)]);
		ieigval	= diag([ones(1,Q)/c1 ones(1,D-Q)/c2]);

		psurf.covars(:,:,m)	= eigvec*eigval*eigvec';
		psurf.Cinv(:,:,m)	= eigvec*ieigval*eigvec';
		psurf.Cidet2(m)		= prod(diag(ieigval));
	end	% m
else
	for m=1:M
		psurf.covars(:,:,m)	= eye(D)/beta(m);
	end
end
%_______________________________________________________________________________
end
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_orient.m in demo mode');

%	create 1st quadrant
t	= 0:.05:pi/2;
x	= [sin(t);cos(t)];
N	= length(x);

%	initialize pseudo psurf
psurf.D	= 2;
psurf.Q	= 1;
psurf.M	= N;
psurf.f	= x;
psurf.x	= t;
psurf.beta	= 10;
psurf.covars	= ones(1,N)/psurf.beta;

%	index of middle node
m		= round(N/2);

a	= [.1 .5 1.5 1.9];
%[row,col]	= lans_fit2page(2*length(a),2);
row=1;col=4;
for sub=1:length(a)
	% orient it
	pstr	= sprintf('-algo 4 -beta 1 -attenuate %0.4f -orient 1',a(sub));
	opsurf	= lans_orient(psurf,pstr);
	% plot mid point only
	subplot(row,col,sub);
	lans_plotmd(opsurf.f,'k-');
	lans_plotmd(opsurf.f(:,m),'ko','-hold 1');
	% unoriented
	lans_plotcovars(psurf.covars(:,m),psurf.f(:,m),'b--',40);
	% oriented
	lans_plotcovars(opsurf.covars(:,:,m),opsurf.f(:,m),'r-',40);
	axis([0 1.3 0 1.3]);
	axis off;
	tstr	= sprintf('\\alpha = %1.2f',a(sub));
	title(tstr);
end
%__________ DEMO ends __________________________________________________________
end