lans_proj2line.m

模式识别工具包

%	lans_proj2line	- Project points onto a piecewise linear curve
%
%	[py,px,result]	= lans_proj2line(y,f,x <,options> <,startx>)
%
%	_____OUTPUTS____________________________________________________________
%	py	projected locations in data space		(col vectors)
%	px	interpolated projections in latent space	(row vector)
%	result	optional results				(cell)	
%		result{1} : squared error		1 x N
%	options							(string)
%		-clos	{0,1}	endpoints are wrapped for closed curves
%			0	default
%			1	connect end points
%	startx	optional index vector for first point if wrapped (col vector)
%
%	_____INPUTS_____________________________________________________________
%	y	test data points in data space			(col vectors)
%	f	ordered points on curve in data space		(col vectors)
%	x	indices of f in latent space			(col vectors)
%
%	_____EXAMPLE____________________________________________________________
%
%
%	_____NOTES______________________________________________________________
%	for demo, call function without parameters
%	- latent variable can be of > 1-D
%
%	_____SEE ALSO___________________________________________________________
%	lans_project
%
%
%	(C) 1999.11.15 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/

function	[py,px,result]	= lans_proj2line(y,f,x,options,startx)
if nargin>0
%__________ REGULAR ____________________________________________________________
if nargin<4
	options	= [];
end
clos	= paraget('-clos',options,0);
val	= paraget('-val',options,0);

[D N]	= size(y);		% N	= # of data points
[D M]	= size(f);		% M	= # of knots

%The following is used instead of ZERO due to difference in precision
%across platforms
%ALSO, in determining smallest distance, if difference between min knot
%and min segment to point is less than prec, we use the seg dist

%prec	= 1e-7;
prec	= 0;

%---------- compute length of closing segment
if clos~=0
	M2	= M+1;
	f(:,M2)	= f(:,1);
	if nargin==5
		x(:,M2)	= startx;
	else
		x(:,M2)	= x(:,1);
	end
else
	M2	= M;
end

%---------- project points onto nearest interval to form M,nf
for p	= 1:N
	y1	= y(:,p);
	ry1	= y1*ones(1,M2);		% replicate x M times

	%-----	Compute dist to each knot (in order)
	f2y1	= ry1-f;			% vector:  knots to point 1:M
	kdist2	= vdist2(f2y1);			% sqr dist to all nodes	 1:M		
	%-----	Find closest KNOT!!!
	minkdist2= min(kdist2(1:M));		% smallest dist to nodes

	%-----	Compute unit segment vector
	l2r	= f(:,2:M2)-f(:,1:M2-1);	% compute segment vector
	dl2r	= vdist(l2r);			% compute segment length
	if dl2r>0
		l2r1	= l2r./(ones(D,1)*dl2r);% unit segment vector
	else
		l2r1	= l2r;
	end

	%-----	Compute projections onto each segment and find valid ones
	projd	= sum(f2y1(:,1:M2-1).*l2r1);	% projected dist on all segments

	vidx1	= find(projd>prec);			% keep ones on right
	vidx2	= find(projd(vidx1)<dl2r(vidx1));	% keep those on segment
	vidx	= vidx1(vidx2);

	if ~isempty(vidx)
		% nearest lies on a segment
		onseg	= (ones(D,1)*projd(vidx)).*l2r1(:,vidx); % proj. vector
		onsegv	= onseg+f(:,vidx); % pos of proj. vector w.r.t. origin
		dist22x	= vdist2(onsegv,y1*ones(1,length(vidx)));
		minsdist2= min(dist22x);
		minsidx	= min(find(dist22x==minsdist2));
		if minsdist2<minkdist2
			% on segment
			py(:,p)	= onsegv(:,minsidx);
			idx	= vidx(minsidx);
			% interpolate
			lambda	= projd(idx)/dl2r(idx);
			px(:,p)	= (1-lambda)*x(:,idx)+lambda*x(:,idx+1);
		else
			% node is nearer
			minkidx	= min(find(kdist2==minkdist2));
			py(:,p)	= f(:,minkidx);
			px(:,p)	= x(:,minkidx);
		end
	else
		% nearest is the node
		minkidx	= min(find(kdist2==minkdist2));
		py(:,p)	= f(:,minkidx);
		px(:,p)	= x(:,minkidx);
	end
	
end	%p

result{1}	= vdist2(y,py);

%__________ REGULAR ends _______________________________________________________
else
%__________ DEMO _______________________________________________________________
clf;clc;
disp('running lans_proj2line.m in demo mode');

y	= randn(3,10);
f	= lans_sphere(0,10,'-linear 1');
x	= sum(lans_cart2pose(f));

lans_plotmd(f,'k-',y,'b*');


[py,px,result]	= lans_proj2line(y,f,x,'-clos 0');
hopen	= lans_plotproj(py,y,'ro-');
hold on
for i=1:length(px)
	text(py(1,i),py(2,i),py(3,i)+.3,sprintf('%3.0f',px(i)));
end

[py,px,result]	= lans_proj2line(y,f,x,'-clos 1',360);
hclose	= lans_plotproj(py,y,'bo-');
hold on
for i=1:length(px)
	text(py(1,i),py(2,i),py(3,i)+.3,sprintf('%3.0f',px(i)));
end

legend([hopen hclose],'open','closed');
%__________ DEMO ends __________________________________________________________
end