lans_ontriplane.m

模式识别工具包

%	lans_ontriplane	- Check if points project ONTO a triangular plane
%
%	[marked,yproj,coef]= lans_ontriplane(y,tri);
%
%	_____OUTPUTS____________________________________________________________
%	marked	denoting whether a point projects onto triplane	(row vector)
%		1:yes 0:no
%	yproj	points on triangle perpendicular to y(:,marked) (col vectors)
%	coef	(all) coefficients of tri approximating y	(col vectors)
%
%	_____INPUTS_____________________________________________________________
%	y	candidate points				(col vectors)
%	tri	triangle defined by 3 col vectors		(col vectors)
%
%	_____NOTES______________________________________________________________
%	- for demo, call function without parameters
%	- if y projects ONTO tri, then it's projection can be express as a
%	  convex combination of the 3 points of tri
%	- ONTO means that projection line is perpendicular to triangular plane 
%	  and projected point is bounded by the triangle
%
%	_____SEE ALSO___________________________________________________________
%	lans_proj2tri
%
%	(C) 1999.06.07 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	[marked,yproj,coef]= lans_ontriplane(y,tri);

if nargin>0
%__________ REGULAR ____________________________________________________________
[d,N]	= size(y);

v31	= tri(:,1)-tri(:,3);
v32	= tri(:,2)-tri(:,3);
v3	= tri(:,3);

A	= [v31 v32];
yp	= y-v3*ones(1,N);

coef		= pinv(A)*yp;
coef(3,:)	= 1-sum(coef);

marked		= ones(1,N);		% assume all points ABOVE triplane
[dum,noton]	= find(coef<0);		% find those that are not

marked(noton)	= zeros(1,length(noton));

pos			= find(marked==1);
thesum			= sum(coef(:,pos));
toobig			= find(thesum>1);
marked(pos(toobig))	= zeros(size(toobig));

remaining		= find(marked==1);
if ~isempty(remaining)
	yproj		= tri*coef(:,remaining);
else
	yproj		= [];
end

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

D	= 3;			% dimension
tri	= rand(D,3);		% generate random triangle
y	= rand(D,20);		% generate random points

[marked,yproj,coef]	= lans_ontriplane(y,tri);

clf;
%patch(tri(1,:),tri(2,:),tri(3,:),'cyan');
lans_plotmd(tri(:,[2 1 3]),'r-','-hold 1');
lans_plotmd(tri(:,[2 3]),'r:','-hold 1');
yes	= find(marked==1);
no	= find(marked==0);
lans_plotmd(y(:,yes),'b*','-hold 1');
lans_plotmd(y(:,no),'ro','-hold 1');

lans_plotproj(yproj,y(:,yes),'m-');

% compute perpendicular vector
if ~isempty(yes)
	pvec	= y(:,yes)-yproj;
	paravec	= tri(:,1)*ones(1,length(yes))-yproj;
	disp('The dotproduct between plane and perpendicular vector');
	dotprod	= sum(paravec.*pvec)	% approx. 0
	plane_se= vdist2(pvec)
end

% projection onto grid line gives larger se
if ~isempty(yes)
	ty	= y(:,yes);
	psurf.f	= [tri tri(:,1)];
	psurf.x	= lans_1speed(psurf.f);
	[prsurf,dum,mse,grid_se]	= lans_project(ty,psurf,'-rsort 0');
	grid_se
	lans_plotproj(prsurf.f,ty,'b-');
end

view(2)
rotate3d on;
display('Rotate figure with mouse')
%__________ DEMO ends __________________________________________________________
end