polyxpoly.m

麻省理工学院的人工智能工具箱,很珍贵,希望对大家有用!

% figure; hold; plotpoly(x1,y1,'b*-'); plotpoly(x2,y2,'ro--')
% [xi yi ii]
% keyboard

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xi,yi,ii,ri] = filterintpts(x1,y1,x2,y2,xi,yi,ii)
%FILTERINTPTS  Filter polygon intersection points.
%   [XI,YI,II,RI] = FILTERINTPTS(X1,Y1,X2,Y2,XI,YI,II) filters the intersection 
%   points from INTPTSALL for the polygon boolean operations.

%  Written by:  A. Kim


err = eps*1e5;

% extract unique intersection points
[a,i,j] = uniqueerr(flipud([xi yi]),'rows',err);
i = length(xi)-flipud(sort(i))+1;
ixy0 = [ii(i,:) xi(i) yi(i)];

% sort all unique intersection points along polygon 1
i1 = sort(unique(ixy0(:,1)));
for n=1:length(i1)
	indx1 = find(ixy0(:,1)==i1(n));
	dist = sqrt((x1(i1(n))-ixy0(indx1,3)).^2 + ...
					(y1(i1(n))-ixy0(indx1,4)).^2);
	[dist,isort] = sort(dist);  indx2 = indx1(isort);
	ixy0(indx1,:) = ixy0(indx2,:);
end

% identify vertex intersection points
nvtx(:,1) = ismembererr(ixy0(:,3:4),[x1 y1],'rows',err);
nvtx(:,2) = ismembererr(ixy0(:,3:4),[x2 y2],'rows',err);
ivtx = find(nvtx(:,1) | nvtx(:,2));

% procedure for determining which vertex intersection points to add:
%	1.  for each vertex intersection point, determine whether the points 
%		before and after the intersection points on polygon 1 are inside
%		or outside polygon 2
%   2.  check conditions to add point

ri = [];  % initialize reverse intersection points index

if ~isempty(ivtx)

% polygon 1 point inside/outside of polygon 2
% for each vertex intersection point, check the last and next midpoints
% on polygon 1 to see if they are inside of polygon 2
	for n=1:length(ivtx)
% determine last point before current vertex intersection point
% (will either be the last point on polygon 1 or an intersection point)
		ilastpt1 = ixy0(ivtx(n),1);
		lastpt = [x1(ilastpt1) y1(ilastpt1)];	% on polygon 1
		if ivtx(n)>1
			if ixy0(ivtx(n)-1,1)==ilastpt1
				lastpt = ixy0(ivtx(n)-1,3:4);	% intersection point
			end
		end
% determine next point after current vertex intersection point
% (will either be the next point on polygon 1 or an intersection point)
		if nvtx(ivtx(n),1)
			inextpt1 = ixy0(ivtx(n),1) + 2;
		else
			inextpt1 = ixy0(ivtx(n),1) + 1;
		end
		nextpt = [x1(inextpt1) y1(inextpt1)];	% on polygon 1
		if ivtx(n)<length(ixy0(:,1))
			if ixy0(ivtx(n)+1,1)==inextpt1-1
				nextpt = ixy0(ivtx(n)+1,3:4);	% intersection point
			end
		end
% compute midpoints and determine whether inside or outside
		lastmidpt = lastpt + .5*(ixy0(ivtx(n),3:4) - lastpt);
		nextmidpt = nextpt + .5*(ixy0(ivtx(n),3:4) - nextpt);
% evaluate last and next points at 9 points: [nw,n,ne;w,c,e;sw,s,se]
		xa = lastmidpt(:,1);  ya = lastmidpt(:,2);
		xb = nextmidpt(:,1);  yb = nextmidpt(:,2);
		p1in(n,1) = inpolygonerr(xa,ya,x2,y2,err);
		p1in(n,2) = inpolygonerr(xb,yb,x2,y2,err);
% 		xca = lastmidpt(:,1);  yca = lastmidpt(:,2);
% 		xcb = nextmidpt(:,1);  ycb = nextmidpt(:,2);
% 		xa = zeros(3,3);  ya = zeros(3,3);
% 		xb = zeros(3,3);  yb = zeros(3,3);
% 		xa(:,1) = xca - xca*eps;  xa(:,2) = xca;  xa(:,3) = xca + xca*eps;
% 		ya(1,:) = yca + yca*eps;  ya(2,:) = yca;  ya(3,:) = yca - yca*eps;
% 		xb(:,1) = xcb - xcb*eps;  xb(:,2) = xcb;  xb(:,3) = xcb + xcb*eps;
% 		yb(1,:) = ycb + ycb*eps;  yb(2,:) = ycb;  yb(3,:) = ycb - ycb*eps;
% 		ina = inpolygon(xa,ya,x2,y2);
% 		inb = inpolygon(xb,yb,x2,y2);
% 		if all(ina(:))
% 			ina = 1;
% 		elseif ~any(ina(:))
% 			ina = 0;
% 		else
% 			ina = .5;
% 		end
% 		if all(inb(:))
% 			inb = 1;
% 		elseif ~any(inb(:))
% 			inb = 0;
% 		else
% 			inb = .5;
% 		end
% 		xa = lastmidpt(:,1);  ya = lastmidpt(:,2);
% 		xb = nextmidpt(:,1);  yb = nextmidpt(:,2);
% 		p1in(n,:) = [ina inb];
	end

% determine which simple (no borders) vertex intersection points to add
%	- check whether points before and after vertex are same or different
% add point for same
	iadd1 = find(ismember(p1in,[0 0; 1 1],'rows'));
% 	nbrdr1 = zeros(size(iadd1));

% determine which complex (borders) vertex intersection points to add:
% 	- check whether points before and after border are same or different
%	- check whether number of borders are odd or even
% add point: diff & odd, same & even
	iadd2 = [];
	indx = [find(p1in(:,2)==.5) find(p1in(:,1)==.5)];
	if ~isempty(indx)
		if length(indx(:,1))==1
			ichk = indx;
		elseif length(indx(:,1))>1
			itmp = find(diff(indx(:,1))>1);
			if isempty(itmp)
				ichk = [1 length(p1in(:,1))];
			else
				ichk = [1 indx(itmp(1),2)];
				for n=1:length(itmp)-1
					ichk = [ichk; indx(itmp(n)+1,1) indx(itmp(n+1),2)];
				end
				ichk = [ichk; indx(itmp(end)+1,1) indx(end)];
			end
		end
		binout = [p1in(ichk(:,1),1) p1in(ichk(:,2),2)];
		bnum = ichk(:,2) - ichk(:,1);
		ipts = find( (binout(:,1)==~binout(:,2) & mod(bnum,2)) | ...
					  (binout(:,1)==binout(:,2) & ~mod(bnum,2)) );
		if ~isempty(ipts)
			iadd2 = ichk(ipts,2);
		end
	end

% procedure for adding additional vertex intersection points:
%	1.  for each vertex intersection point to be added, determine ordering
%	2.  add points with correct ordering

% determine ordering for non-border vertex intersection points:
% for each non-border vertex intersection point to add, check the last and 
% next midpoints on polygon 2 to see if they are inside of polygon 1 to 
% determine ordering of intersection points
	rvsordr1 = [];
	if ~isempty(iadd1)
		for n=1:length(iadd1)
% specify current vertex point to add
			ivtxpt = ivtx(iadd1(n));
% determine last point before current vertex intersection point or border
% point (last point on polygon 2 or an intersection point)
			ilastpt2 = ixy0(ivtxpt,2);
			lastpt = [x2(ilastpt2) y2(ilastpt2)];	% on polygon 2
			if ivtxpt>1
				if ixy0(ivtxpt-1,2)==ilastpt2
					lastpt = ixy0(ivtxpt-1,3:4);	% intersection point
				end
			end
% determine next point after current vertex intersection point
% (will either be the next point on polygon 1 or an intersection point)
			if nvtx(ivtxpt,2)
				inextpt2 = ixy0(ivtxpt,2) + 2;
			else
				inextpt2 = ixy0(ivtxpt,2) + 1;
			end
			nextpt = [x2(inextpt2) y2(inextpt2)];	% on polygon 2
			if ivtxpt<length(ixy0(:,1))
				if ixy0(ivtxpt+1,2)==inextpt2-1
					nextpt = ixy0(ivtxpt+1,3:4);	% intersection point
				end
			end
% compute midpoints and determine whether inside or outside
			lastmidpt = lastpt + (ixy0(ivtxpt,3:4) - lastpt) / 2;
			nextmidpt = nextpt + (ixy0(ivtxpt,3:4) - nextpt) / 2;
			xa = lastmidpt(:,1);  ya = lastmidpt(:,2);
			xb = nextmidpt(:,1);  yb = nextmidpt(:,2);
% 			p2in(n,:) = inpolygon([xa xb],[ya yb],x1,y1);
			p2in(n,1) = inpolygonerr(xa,ya,x1,y1,err);
			p2in(n,2) = inpolygonerr(xb,yb,x1,y1,err);
		end
% if polygon 2 is counter-clockwise (intersection and union), reverse 
% order if similar in/out states; if polygon 2 is clockwise (subtraction),
% reverse oder if opposite in/out states
		if sparea(x2,y2)<0
			rvsordr1 = p1in(iadd1,1)==p2in(:,1);
		else
			rvsordr1 = p1in(iadd1,1)~=p2in(:,1);
		end
	end

% determine ordering for border vertex intersection points:
%	- if the point on polygon 2 before the current vertex intersection point 
%	  is on the border of polygon 1, then the two polygons follow the same 
% 	  direction along their common border
%	- if the point on polygon 2 before the current vertex intersection point  
%	  is not on the border of polygon 1, then the two polygons do not follow  
%	  the same direction along their common border
	rvsordr2 = [];
	if ~isempty(iadd2)
		for n=1:length(iadd2)
% specify current vertex point to add
			ivtxpt = ivtx(iadd2(n));
% point on polygon 2 before the current vertex intersection point
			ipt = ixy0(ivtxpt,1:2);
% reverse order border segments are in opposite direction
			dir1 = atan2(y1(ipt(1)+1)-y1(ipt(1)),x1(ipt(1)+1)-x1(ipt(1)));
			dir2 = atan2(y2(ipt(2)+1)-y2(ipt(2)),x2(ipt(2)+1)-x2(ipt(2)));
% 			if dir1~=dir2
			if abs(dir1-dir2)>err
				rvsordr2(n,1) = 1;
			else
				rvsordr2(n,1) = 0;
			end
		end
% 		rvsordr2 = inpolygon(x2(ipt),y2(ipt),x1,y1)~=.5;
	end

% add points using add points index and reverse order flag
	[iaddpts,isort] = sort(ivtx([iadd1; iadd2]));
	rvsordr = [rvsordr1; rvsordr2];
	rvsordr = rvsordr(isort);
	ixy = ixy0;
	for n=1:length(iaddpts)
		iins = find(ismembererr(ixy(:,3:4),ixy0(iaddpts(n),3:4),'rows',err));
		iadd = find(ismembererr([xi yi],ixy0(iaddpts(n),3:4),'rows',err));
		if rvsordr(n)
			if length(iadd(:)) > 1
				iadd = flipud(iadd(1:2));
			end
			ri = [ri; iins];  % store reverses
		else
			iadd = iadd(1:2);
		end
		addvtxpts = [ii(iadd,:) xi(iadd) yi(iadd)];
		ixy = [ixy(1:iins-1,:); addvtxpts; ixy(iins+1:end,:)];
	end

% all intersection points
	ixy0 = ixy;	

end

% set outputs
ii = ixy0(:,1:2);  xi = ixy0(:,3);  yi = ixy0(:,4);

% check for epsilon error
[xi,yi] = ptserr(xi,yi,[x1;x2],[y1;y2],err);