lans_proj2sqr.m

模式识别工具包

%	lans_proj2sqr	- Project data onto a square patch in high-D space
%
%	[se,np,idx,gse]	= lans_proj2sqr(y,sqrv,sqri[,option])
%
%	_____OUTPUTS____________________________________________________________
%	se	squared error from each vector y to patch	(row vector)
%	np	nearest point on manifold			(col vectors)
%	idx	index of np w.r.t. sqri				(col vectors)
%	gse	Grid based squared error			(row vector)
%
%	_____INPUTS_____________________________________________________________
%	y	data						(col vectors)
%	sqrv	4 vector points of the square patch ordered as	(col vectors)
%			1 2
%			3 4
%		
%	sqri	manifold indices of the 4 points		(col vectors)
%	option	interpolation type				(string)	
%		'grid'
%		'plane'
%		'triangle'
%
%	_____NOTES______________________________________________________________
%	- for printable demo, call function without parameters
%	- only works for 2-D sub-manifolds
%	- se is squared error, NOT distance!
%	- Euclidean distance computed from sqrv to interpolate sqri
%	- 'grid'
%	  only the principal 'right' and 'down' grids will be considered, an
%	  external pass is needed to take care of the uncovered external
%	  frame areas
%	- 'plane'	(NOT RECOMMENDED for INTERPOLATION)
%	  assumes principal 'right' 12 and 'down'13 grid vectors are approx.
%	  perpendicular and all 4 points approx. lie on a plane.
%	  Interpolation based on right and down grid vectors.
%	- 'triangle'
%	  both triangulations are tried and the min result
%	  recorded. In addition, a 'grid' pass is run because the projection
%	  point obtained via 'grid' may be more accurate than that obtained
%	  via solving the pseudo-inverse of the triangle
%
%	- gse used in expermental evaluation, only works with 'triangle' option
%
%	_____SEE ALSO___________________________________________________________
%	lans_proj2tri lans_psurfproj
%
%	(C) 1999.11.10 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______________________________________________________________

function	[se,np,idx,gse]	= lans_proj2sqr(y,sqrv,sqri,option,debug)

if nargin>0
%__________ REGULAR ____________________________________________________________

if nargin<5
	debug=0;
	if nargin<4
		option	= 'plane';
	end
end
gse	= 0;

[D,N]	= size(y);
[Q,dum]	= size(sqrv);

switch(option)
case 'grid'	
	pvec	= y-sqrv(:,1)*ones(1,N);		% point vector
	plen	= vdist(pvec);

	%______	Initialize nearest projection to patch origin node
	se	= plen;					% dist to nearest node
	np	= zeros(D,N);				% origin shifted to pvec
	idx	= sqri(:,1)*ones(1,N);			% nearest node idx

	% compute vectors
	rvec	= sqrv(:,2)-sqrv(:,1);			% right vector
	rlen	= norm(rvec);
	if rlen~=0
		rvec1	= rvec/rlen;			% unit right vector
	else
		rvec1	= rvec;
	end

	dvec	= sqrv(:,3)-sqrv(:,1);			% down vector
	dlen	= norm(dvec);
	if dlen~=0
		dvec1	= dvec/dlen;			% unit down vector
	else
		dvec1	= dvec;
	end

	%_____	Projection onto right vector
	rproj	= rvec1'*pvec;				% right projected len
	m1	= find(rproj>=0 & rproj<rlen);		% within rvec
	if ~isempty(m1)
		np(:,m1)= rvec1*rproj(m1);		% valid proj right_vec
		if rlen~=0
			a	= rproj(m1)/rlen;
		else
			a	= rproj(m1);
		end
		% assign interpolated grid points
		se(m1)		= plen(m1).*plen(m1)-rproj(m1).*rproj(m1);
		idx(:,m1)	= sqri(:,1)*(1-a) + sqri(:,2)*a;
	end

	%_____	Projection onto down vector
	dproj	= dvec1'*pvec;				% down projected len
	m2	= find(dproj>=0 & dproj<dlen);		% within dvec
	if ~isempty(m2)
		p2dvec	= dvec1*dproj(m2);		% valid proj down_vec
		% assign interpolated grid points if smaller than ontoright
		dse	= plen(m2).*plen(m2)-dproj(m2).*dproj(m2);
		wsmall	= find(dse<se(m2));
	
		if ~isempty(wsmall)
			gidx		= m2(wsmall);		% global index
			if dlen~=0
				a		= dproj(gidx)/dlen;
			else
				a		= dproj(gidx);
			end
			se(gidx)	= dse(wsmall);
			np(:,gidx)	= p2dvec(:,wsmall);
			idx(:,gidx)	= sqri(:,1)*(1-a) + sqri(:,3)*a;
		end
	end
	np	= np + sqrv(:,1)*ones(1,N);		% shift back to origin

case 'plane'	
	% assumes right and down vector approximately perpendicular
	pvec	= y-sqrv(:,1)*ones(1,N);		% point vector
	plen	= vdist(pvec);

	rvec	= sqrv(:,2)-sqrv(:,1);			% right vector
	rlen	= norm(rvec);
	if rlen~=0
		rvec1	= rvec/rlen;			% unit right vector
	else
		rvec1	= rvec;
	end

	dvec	= sqrv(:,3)-sqrv(:,1);			% down vector
	dlen	= norm(dvec);
	if dlen~=0
		dvec1	= dvec/dlen;			% unit down vector
	else
		dvec1	= dvec;
	end

	%______	Initialize nearest projection to patch origin node
	se	= plen;					% dist to nearest node
	np	= zeros(D,N);				% origin shifted to pvec
	idx	= sqri(:,1)*ones(1,N);			% nearest node idx

	% project
	rproj	= rvec1'*pvec;				% right projections
	win	= find(rproj>=0 & rproj<=rlen);		% winner index (x dir)
	
	pvec2	= pvec(:,win);
	plen2	= plen(:,win);

	dproj	= dvec1'*pvec2;				% down projections (sub)
	win2	= find(dproj>=0 & dproj<=dlen);		% winner index (y dir)
	gwin	= win(win2);				% winner index (global)
	
	if ~isempty(gwin)	
		np(:,gwin)	= rvec1*rproj(gwin) + dvec1*dproj(win2);
		nplen		= vdist2(np(:,gwin));

		% interpolated indices
		if rlen~=0
			a		= rproj(gwin)/rlen;
		else
			a		= rproj(gwin);
		end
		if dlen~=0
			b		= dproj(win2)/dlen;
		else
			b		= dproj(win2);
		end

		idx(:,gwin)	= sqri(:,1)*(1-a) + sqri(:,2)*a + sqri(:,1)*(1-b)+sqri(:,3)*b;
		se(gwin)	= plen(gwin).*plen(gwin) - nplen;
	end

	np	= np + sqrv(:,1)*ones(1,N);		% shift back to origin

case 'triangle'	
%	4 triangulations:
%		11 2	3 44
%		1 22	33 4

	%______	Initialize nearest projection to patch origin node
	np	= sqrv(:,1)*ones(1,N);			% nearest = origin node
	idx	= sqri(:,1)*ones(1,N);			% nearest node idx
	pvec	= y-np;				% point vector
	nnse	= vdist(pvec);			% nearest node distance
	se	= nnse;

	%_____	The 4 triangulations
	triv	= {sqrv(:,[1 2 3]),sqrv(:,[4 3 2]),sqrv(:,[2 1 4]),sqrv(:,[3 4 1])};
	tridx	= {sqri(:,[1 2 3]),sqri(:,[4 3 2]),sqri(:,[2 1 4]),sqri(:,[3 4 1])};

	% for each triangle, proj y onto it
	for t=1:4
		thetri		= triv{t};
		thetridx	= tridx{t};
		[tse,tnp,tidx]	= lans_proj2tri(y,thetri,thetridx);
		ws		= find(tse<se);
		if ~isempty(ws)
			se(ws)		= tse(ws);
			np(:,ws)	= tnp(:,ws);
			idx(:,ws)	= tidx(:,ws);
		end
	end

	% project all y onto grid in case it is nearer
	% not possible in theory, but in practice numerical errors associated
	% with solving the plane projection coefficients result in the grid
	% projection giving lower errors
 	[gse,snp,sidx]	= lans_proj2sqr(y,sqrv,sqri,'grid',debug);
	ws		= find(gse<se);
	if ~isempty(ws)
		se(ws)		= gse(ws);
		np(:,ws)	= snp(:,ws);
		idx(:,ws)	= sidx(:,ws);
	end
end	% case
%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_proj2sqr.m in demo mode');

if ~exist('rseed')
	rseed	= now;
end
rand('state',rseed);

N	= 1;
D	= 3;
y	= [.2 .7 .6]';

sv	= [0 1 0;1 1 1;0 0 0;1 0 0]';
si	= [1 2 3 4];

ptype={'nn','grid','triangle1'};
plabel={'MSE_{nn}','MSE_{grid}','MSE_{\Delta}','MSE_{\Delta2}'};
number={'a','b','c','d'};
figure(1)
clf
tf	= length(ptype);

for f=1:tf
	subplot(1,3,f);
	if f==4
		[se,np,idx]	= lans_proj2tri(y,sv(:,[1 2 4]),si([1 2 4]));	
	end
	if f==3
		[se,np,idx]	= lans_proj2tri(y,sv(:,[1 2 3]),si([1 2 3]));	
	end
	if f==2
		[se,np,idx]	= lans_proj2sqr(y,sv,si,char(ptype{f}));
	end
	if f==1
		nndist	= vdist2(y*ones(1,4),sv);
		se	= min(nndist);
		idx	= find(nndist==se);
		np	= sv(:,idx);
	end

	lans_plotmd(sv(:,[2 1]),'ro-','-hold 1');
	lans_plotmd(sv(:,[3 1]),'ro-','-hold 1');
	lans_plotmd(sv(:,[2 4 3]),'ro--','-hold 1');
	lans_plotmd(y,'b*','-hold 1');
	lans_plotproj(np,y,'k--');
	rotate3d on;
	for i=1:size(y,2)
		text(y(1,i),y(2,i),y(3,i)+.3,'\bf{y}');	
	end
	mse	= mean(se);
	tstr	= sprintf('(%s) %s = %0.4f',char(number{f}),char(plabel{f}),mse);
	title(tstr);

	if (f==3)
		fill3(sv(1,1:3),sv(2,1:3),sv(3,1:3),'c')
	end
	if (f==4)
		fill3(sv(1,[1 2 4]),sv(2,[1 2 4]),sv(3,[1 2 4]),'c')
	end
	axis off;
	view(15,30)
end
%__________ DEMO ends __________________________________________________________
end