lans_psurfproj.m

模式识别工具包

%	lans_psurfproj	- Project data (approx.) onto principal surface
%
%	[prsurf,mse,se,gse]	= lans_psurfproj(y,psurf[,itype])
%
%	_____OUTPUTS____________________________________________________________
%	prsurf	subset of psurf					(psurf)
%		.idx	indices w.r.t. psurf if uninterpolated	(row vector)
%		.x	interpolated projected knots		(col vectors)
%		.f	interpolated projected manifold values	(col vectors)
%	xs	estimated projection indices			(row vector)
%		e.g. psurf.x(pidx) gives the projected knots
%		or
%		interpolated indices
%	mse	mean squared error of data to psurf.f		(scalar)
%	se	squared error/dist of each data to psurf.f	(row vector)
%	gse	squared error for 'grid' used for evaluation	(row vector)
%		works only when itype='plane'
%
%	_____INPUTS_____________________________________________________________
%	y	unormalized data in original dimension		(col vectors)
%	psurf	Q-D principal surface				(psurf)
%		see lans_psurfinit
%	itype	projection interpolation type for >1-D manifold	(string)
%		nn	distance to nearest neighboring knot	(default)
%		grid	distance to nearest grid line (points inclusive)
%		plane	distance to nearest plane
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters
%	- for 'plane' type, all 4 triangulations of a patch are tried to find
%	  the nearest projection
%	- if point does not project directly onto manifold surface,
%	  nearest grid projection is located.
%	- SORTING is not applicable for >1-D manifold
%	- define psurf.nn to use nearest neighbor projection for manifold
%	dimension Q>1, where nearest latent vector is designated as
%	projection point, otherwise stringwise approximation is used.
%	- prsurf.idx is not applicable with interplation (e.g. Q=1)
%
%	_____SEE ALSO___________________________________________________________
%	lans_project	lans_proj2sqr
%
%	(C) 1999.05.28 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______________________________________________________________
%	- let 'plane' also output 'grid' error (for experiment evaluation)

function	[prsurf,mse,se,gse]	= lans_psurfproj(y,psurf,itype)

if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin<3
	itype	= 'nn';
end

gse	= 0;

[D N]	= size(y);

M	= psurf.M;
Q	= psurf.Q;

%__________	1-D manifolds (i.e. line)
if Q==1
	if strcmp(itype,'nn')
		distM		= lans_mahadist(psurf.f,y);
		[sdistM,idx]	= sort(distM);
		rse		= sdistM(1,:);
		se		= rse.*rse;
		prsurf.idx	= idx(1,:);
		prsurf.x	= psurf.x(:,idx(1,:));
		prsurf.f	= psurf.f(:,idx(1,:));
		mse		= mean(se);
	else
		[prgsurf,prgdata,mse,se]	= lans_project(y,psurf,'-rsort 0');
		prsurf.idx	= [];
		prsurf.x	= prgsurf.x;
		prsurf.f	= prgsurf.f;
	end
%__________	2-D manifolds
elseif Q==2
    M1	= round((psurf.M)^(1/psurf.Q));		% # nodes in each dimension
    switch(itype)	
    %__________	Simple nearest neighbor is used for >1-D Manifolds
    case 'nn',
	distM		= lans_mahadist(psurf.f,y);
	[sdistM,idx]	= sort(distM);
	rse		= sdistM(1,:);
	se		= rse.*rse;
	prsurf.idx	= idx(1,:);
	prsurf.x	= psurf.x(prsurf.idx);
	prsurf.f	= psurf.f(:,prsurf.idx);
    %__________	compute projection distance to nearest gridline
    case 'grid',
	xidx		= psurf.xidx;		% indices to 2-D manifold
	%_____ project all points onto bottom-right-most grid not in patch grids
	% i.e. mirror image of 'L'
	invLidx		= [xidx(end,1:end) xidx(end-1:-1:1,end)'];
	bsurf.f		= psurf.f(:,invLidx);
	bsurf.x		= lans_1speed(bsurf.f);
	bsurf.xQ	= psurf.x(:,invLidx);
	[bprsurf,dum,dum,minse]	= lans_project(y,bsurf,'-rsort 0');
	minnp		= bprsurf.f;
	minidx		= bprsurf.xQ;

	%_____	project all points onto all (M1-1)^2 patch grids
	for x1=1:M1-1
		for x2=1:M1-1
			patchidx= xidx([x1 x1+1],[x2 x2+1])';
			sv	= psurf.f(:,patchidx);
			si	= psurf.x(:,patchidx);
			[se,np,idx]= lans_proj2sqr(y,sv,si,'grid');	
			ws	= find(se<minse);
			minnp(:,ws)	= np(:,ws);
			minse(ws)	= se(ws);
			minidx(:,ws)	= idx(:,ws);
		end
	end
	prsurf.f	= minnp;
	prsurf.x	= minidx;
	se		= minse;
	mse		= mean(se);
	gse		= se;

    case 'plane'
    	f	= psurf.f;		% points on 2-D manifold
	xidx	= psurf.xidx;		% indices to 2-D manifold

	%_____ project all points onto bottom-right-most grid not on patch grids
	% i.e. mirror image of 'L'
	invLidx		= [xidx(end,1:end) xidx(end-1:-1:1,end)'];
	bsurf.f		= psurf.f(:,invLidx);
	bsurf.x		= lans_1speed(bsurf.f);
	bsurf.xQ	= psurf.x(:,invLidx);
	[bprsurf,dum,dum,minse]	= lans_project(y,bsurf,'-rsort 0');
	minnp		= bprsurf.f;
	minidx		= bprsurf.xQ;

	gse		= minse;
	%_____	project all points onto all (M1-1)^2 patches
	for x1=1:M1-1
		for x2=1:M1-1
			patchidx	= xidx([x1 x1+1],[x2 x2+1])';
			sv		= psurf.f(:,patchidx);
			si		= psurf.x(:,patchidx);
			[se,np,idx,thegse]= lans_proj2sqr(y,sv,si,'triangle');	
			ws		= find(se<minse);

			minse(ws)	= se(ws);
			minnp(:,ws)	= np(:,ws);
			minidx(:,ws)	= idx(:,ws);

			gws		= find(thegse<gse);
			gse(gws)	= thegse(gws);
		end
	end

	prsurf.f	= minnp;
	prsurf.x	= minidx;
	se		= minse;
	mse		= mean(se);

    end		% switch
    mse	= mean(se);
end
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_psurfproj.m in demo mode');

dstring	= '-algo 4 -attenuate 1 -beta 0 -covartype spherical -L 4 -manifold 2square -s 2 -M 25 -orient 0 ';

rseed	= 1;
%_____	Data
N	= 10;
D	= 3;
ball.var	= .1;
ball.radiusOut	= 1.5;
ball.radiusIn	= .5;
ball.center	= zeros(D,1);
y		= lans_genball(N,rseed,ball);
[psurf,prsurf,mse,se]	= lans_psurfinit(y,dstring);

subplot(2,1,1);
for k=1:5
	psurf	= lans_psurf1(y,psurf,dstring);
	lans_plotmd(y,'yo');
	lans_psurfpgrid(psurf,'g-');
	drawnow
	ptype	= {'nn','grid','plane';'g-','b-','r-'};
	pmse	= zeros(1,length(ptype));
	for i=1:3
		flops(0);
		[prsurf,mse,se,gse] = lans_psurfproj(y,psurf,char(ptype{1,i}));
		pflops(i)	= flops;
		lans_plotproj(prsurf.f,y,char(ptype{2,i}));
		pmse(i)		= mse;
	end
	tstr=sprintf('[%d] nn=%0.4f grid=%0.4f plane=%0.4f',k, pmse(1),pmse(2),pmse(3));
	title(tstr);
	drawnow;
	pmse
end

pflops
for n=1:N
	text(y(1,n),y(2,n),y(3,n),num2str(n));
end

rotate3d on;
lans_plotmd(prsurf.f,'m*','-hold 1');

%disp('actual squared distance');
%vdist2(prsurf.f-y)

% grid
%
subplot(2,1,2);
gridsurf	= psurf;
gridsurf.f	= psurf.x;
lans_psurfpgrid(gridsurf);
%lans_plotmd(prsurf.x,'m*','-hold 1');
for n=1:N
	text(prsurf.x(1,n),prsurf.x(2,n),num2str(n));
end

disp('Plane psurfproj shown');
%__________ DEMO ends __________________________________________________________
end