frontpts.m

为了下载东西

function [xf,yf,zf] = frontpts(x,y,z,d,ax,page)
% FRONTPTS  Move points toward (or away from) the viewer.
%
%   This is useful because points or lines plotted exactly on a surface
%   often don't show up (they are seen to be infinitesimally behind the
%   surface), so FRONTPTS can return new points such that they appear
%   to be in the same place but are really just closer to the viewer.
%   FRONTPTS DEMO will display a demo of this.
%
%   [XF,YF,ZF] = FRONTPTS(X,Y,Z,D) moves the (X,Y,Z) points toward the
%   user by a distance D and returns the new points in 3-D space.
%   X, Y, Z, and D can be scalars, vectors or matrices.  The number of
%   rows in each must be the same or unity (in which case it will be
%   replicated as necessary to match the size of the other variables).
%   Likewise with the number of columns.
%
%   Note that D is not the distance in Euclidean space, but rather some
%   measure of how much the points should be moved toward the viewer in
%   a transformed space.
%
%   The axes can be specified with FRONTPTS(...,AX) where AX is a handle
%   to the desired axes.  It defaults to the current axes.
%
%   Spacing on the printed page can be used with FRONTPTS(...,'page').

% Copyright (c)1995, Erik A. Johnson <[e-mail]>, 11/13/95

% Revision history:
%   11/13/95  EAJ  Corrected handling for 'reverse' axis directions
%   10/06/95  EAJ  Corrected occasional conflict with SUBPLOT
%   24.10.97  HK   Color of contour lines set in contour() instead of with set()
%                  Version specific: Matlab 5 uses DataAspectRatio -> AspectRatio                 

% demo
dodemo=0;
if (nargin==0), dodemo=1; elseif ((nargin==1)&(isstr(x))), dodemo=1; end;
if (dodemo),
	% create some data
	[x,y,z] = peaks(21);
	a=.05;
	% set the axis limits
	clf;
	hold on;
	view(viewmtx(-37.5,30,20));
	minx=min([x(:);y(:)]);
	maxx=max([x(:);y(:)]);
	axis([minx maxx minx maxx minx maxx]);
	xlabel('x'); ylabel('y'); zlabel('z'); title('FRONTPTS demo');
	% plot the surface
	hs1=surf(x,ones(size(x))*(minx*(1-a)+maxx*a),y,z);
	hs2=surf(x,ones(size(x))*(minx*a+maxx*(1-a)),y,z);
	set([hs1;hs2],'EdgeColor','none','FaceColor','interp');
	% create and plot the contour lines
	[C,H1] = contour(x,y,z,11,'k');
	[C,H2] = contour(x,y,z,11,'k');
	% move them to the proper place
	for k=1:length(H1),
		xdat = get(H1(k),'XData');
		ydat = ones(size(xdat));
		zdat = get(H1(k),'YData');
		set(H1(k),'XData',xdat,'YData',ydat*(minx*a+maxx*(1-a)),'ZData',zdat);
		ydat=ydat*(minx*(1-a)+maxx*a);
		[xdat,ydat,zdat] = frontpts(xdat,ydat,zdat,1);
		set(H2(k),'XData',xdat,'YData',ydat,'ZData',zdat);
	end;
	% text
	t1=text(maxx*1.03-minx*.03,minx*a+maxx*(1-a),maxx,'Contour lines in plane');
	t2=text(maxx,minx*(1-a)+maxx*a,maxx*1.03-minx*.03,'Contour lines in front of plane','HorizontalAlignment','right');
	hold off;
	return;
end;

if (nargin<4),
	error('FRONTPTS requires at least 4 input arguments.');
elseif (nargout>3),
	error('FRONTPTS produces at most 3 output arguments.');
elseif (nargin>6),
	error('FRONTPTS takes at most 6 input arguments.');
elseif (nargin==4),
	ax = [];
	page = [];
elseif (nargin==5),
	if (isstr(ax)),
		page = 1;
		ax = [];
	else,
		page = [];
	end;
else, %nargin==6
	if (isstr(ax)),
		ax = page;
		page = 1;
	end;
end;

% Set defaults
if (~all(size(ax))), ax=gca; end;
page = all(size(page));

% get number of points and correct for scalar/vector arguments
sizes = [size(x);size(y);size(z);size(d);size(ax);size(page)];
maxsize = max(sizes);
if (any((sizes(:)~=kron(maxsize',ones(6,1)))&(sizes(:)~=1))),
	error('FRONTPTS requires that all inputs be the same size or scalar');
end;
if (maxsize(1)>1),
	o = ones(maxsize(1),1);
	if (size(x,1)==1), x = o*x; end;
	if (size(y,1)==1), y = o*y; end;
	if (size(z,1)==1), z = o*z; end;
	if (size(d,1)==1), d = o*d; end;
	if (size(ax,1)==1), ax = o*ax; end;
	if (size(page,1)==1), page = o*page; end;
end;
if (maxsize(2)>1),
	o = ones(1,maxsize(2));
	if (size(x,2)==1), x = x*o; end;
	if (size(y,2)==1), y = y*o; end;
	if (size(z,2)==1), z = z*o; end;
	if (size(d,2)==1), d = d*o; end;
	if (size(ax,2)==1), ax = ax*o; end;
	if (size(page,2)==1), page = page*o; end;
end;
xf=x; yf=y; zf=z;
pts=[x(:) y(:) z(:)]'; d=d(:)'; ax=ax(:)'; page=page(:)';
narrows = size(pts,2);

% Get axes limits, range, min; correct for aspect ratio and log scale
axm      = zeros(3,narrows);
axr      = zeros(3,narrows);
axrev    = zeros(3,narrows);
ap       = zeros(2,narrows);
xyzlog   = zeros(3,narrows);
limmin   = zeros(2,narrows);
limrange = zeros(2,narrows);
oneax = all(ax==ax(1));
if (oneax),
	T = zeros(4,4);
	invT = zeros(4,4);
else,
	T        = zeros(16,narrows);
	invT     = zeros(16,narrows);
end;
axnotdone = ones(size(ax));
while (any(axnotdone)),
	ii = min(find(axnotdone));
	curax = ax(ii);
	curpage = page(ii);
	% get axes limits and aspect ratio
	axl = [get(curax,'XLim'); get(curax,'YLim'); get(curax,'ZLim')];

   %Allow for different treatment in Matlab 4 and Matlab 5
	vers=version;vers=str2num(vers(1));
	if (vers>4), ar = get(curax,'DataAspectRatio');
   else, ar = get(curax,'AspectRatio');
   end;
	% get axes size in pixels (points)
	u = get(curax,'Units');
	if (curpage),
		curfig = get(curax,'Parent');
		pu = get(curfig,'PaperUnits');
		set(curfig,'PaperUnits','points');
		pp = get(curfig,'PaperPosition');
		set(curfig,'PaperUnits',pu);
		set(curax,'Units','normalized');
		curap = get(curax,'Position');
		curap = pp.*curap;
	else,
		set(ax,'Units','pixels');
		curap = get(curax,'Position');
	end;
	set(curax,'Units',u);
	% adjust limits for log scale on axes
	curxyzlog = [strcmp(get(curax,'XScale'),'log'); ...
	             strcmp(get(curax,'YScale'),'log'); ...
	             strcmp(get(curax,'ZScale'),'log')];
	if (any(curxyzlog)),
		ii = find([curxyzlog;curxyzlog]);
		if (any(axl(ii)<=0)),
			error('FRONTPTS does not support non-positive limits on log-scaled axes.');
		else,
			axl(ii) = log10(axl(ii));
		end;
	end;
	% correct for 'reverse' direction on axes;
	curreverse = [strcmp(get(curax,'XDir'),'reverse'); ...
	              strcmp(get(curax,'YDir'),'reverse'); ...
	              strcmp(get(curax,'ZDir'),'reverse')];
	ii = find(curreverse);
	if (all(size(ii))),
		axl(ii,[1 2])=-axl(ii,[2 1]);
	end;
	% correct for aspect ratio
	if (~isnan(ar(1))),
		if (curap(3) < ar(1)*curap(4)),
			curap(2) = curap(2) + (curap(4)-curap(3)/ar(1))/2;
			curap(4) = curap(3)/ar(1);
		else,
			curap(1) = curap(1) + (curap(3)-curap(4)*ar(1))/2;
			curap(3) = curap(4)*ar(1);
		end;
	end;
	% correct for 'equal'
	% may only want to do this for 2-D views, but seems right for 3-D also
	if (~isnan(ar(2))),
		if ((curap(3)/(axl(1,2)-axl(1,1)))/(curap(4)/(axl(2,2)-axl(2,1)))>ar(2)),
			incr = curap(3)*(axl(2,2)-axl(2,1))/(curap(4)*ar(2)) - (axl(1,2)-axl(1,1));
			axl(1,:) = axl(1,:) + incr/2*[-1 1];
		else,
			incr = ar(2)*(axl(1,2)-axl(1,1))*curap(4)/curap(3) - (axl(2,2)-axl(2,1));
			axl(2,:) = axl(2,:) + incr/2*[-1 1];
		end;
	end;
	% compute the range of 2-D values
	curT = get(curax,'Xform');
	lim = curT*[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1;0 0 0 0 1 1 1 1;1 1 1 1 1 1 1 1];
	lim = lim(1:2,:)./([1;1]*lim(4,:));
	curlimmin = min(lim')';
	curlimrange = max(lim')' - curlimmin;
	curinvT = inv(curT);
	if (~oneax),
		curT = curT.';
		curinvT = curinvT.';
		curT = curT(:);
		curinvT = curinvT(:);
	end;
	% check which arrows to which cur corresponds
	ii = find((ax==curax)&(page==curpage));
	oo = ones(1,length(ii));
	axr(:,ii)      = diff(axl')' * oo;
	axm(:,ii)      = axl(:,1)    * oo;
	axrev(:,ii)    = curreverse  * oo;
	ap(:,ii)       = curap(3:4)' * oo;
	xyzlog(:,ii)   = curxyzlog   * oo;
	limmin(:,ii)   = curlimmin   * oo;
	limrange(:,ii) = curlimrange * oo;
	if (oneax),
		T    = curT;
		invT = curinvT;
	else,
		T(:,ii)    = curT    * oo;
		invT(:,ii) = curinvT * oo;
	end;
	axnotdone(ii) = zeros(1,length(ii));
end;

% correct for log scales
ii = find(xyzlog(:));
if (all(size(ii))),
	pts(ii) = real(log10(pts(ii)));
end;

% correct for reverse directions
jj = find(axrev(:));
if (all(size(jj))),
	pts(jj) = -pts(jj);
end;

% convert to 2-space
tmp1=[(pts-axm)./axr; ones(1,narrows)];
if (oneax), X=T*tmp1;
else, tmp1=[tmp1;tmp1;tmp1;tmp1]; tmp1=T.*tmp1;
      tmp2=zeros(4,4*narrows); tmp2(:)=tmp1(:);
      X=zeros(4,narrows); X(:)=sum(tmp2)'; end;
X = X ./ (ones(4,1)*X(4,:));

% add the distance
X(3,:) = X(3,:) + d;

% convert back to 3-space
if (oneax), threeD=invT*X;
else, tmp1=[X;X;X;X]; tmp1=invT.*tmp1;
      tmp2=zeros(4,4*narrows); tmp2(:)=tmp1(:);
      threeD=zeros(4,narrows); threeD(:)=sum(tmp2)'; end;
pts = (threeD(1:3,:)./(ones(3,1)*threeD(4,:))).*axr+axm;

% recorrect for reverse directions
if (all(size(jj))),
	pts(jj) = -pts(jj);
end;

% convert back to log space
if (all(size(ii))),
	pts(ii) = 10.^(pts(ii));
end;

% put in output matrices
xf(:) = pts(1,:).';
yf(:) = pts(2,:).';
zf(:) = pts(3,:).';