www.gusucode.com > 冈萨雷斯数字图像处理matlab版源码V1.1 > 冈萨雷斯数字图像处理matlab版源码V1.1.3/code/新建文件夹/polyangles.m
function angles = polyangles(x, y) %POLYANGLES Computes internal polygon angles. % ANGLES = POLYANGLES(X, Y) computes the interior angles (in % degrees) of an arbitrary polygon whose vertices are given in % [X, Y], ordered in a clockwise manner. The program eliminates % duplicate adjacent rows in [X Y], except that the first row may % equal the last, so that the polygon is closed. % Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins % Digital Image Processing Using MATLAB, Prentice-Hall, 2004 % $Revision: 1.6 $ $Date: 2003/11/21 14:44:06 $ % Preliminaries. [x y] = dupgone(x, y); % Eliminate duplicate vertices. xy = [x(:) y(:)]; if isempty(xy) % No vertices! angles = zeros(0, 1); return; end if size(xy, 1) == 1 | ~isequal(xy(1, :), xy(end, :)) % Close the polygon xy(end + 1, :) = xy(1, :); end % Precompute some quantities. d = diff(xy, 1); v1 = -d(1:end, :); v2 = [d(2:end, :); d(1, :)]; v1_dot_v2 = sum(v1 .* v2, 2); mag_v1 = sqrt(sum(v1.^2, 2)); mag_v2 = sqrt(sum(v2.^2, 2)); % Protect against nearly duplicate vertices; output angle will be 90 % degrees for such cases. The "real" further protects against % possible small imaginary angle components in those cases. mag_v1(~mag_v1) = eps; mag_v2(~mag_v2) = eps; angles = real(acos(v1_dot_v2 ./ mag_v1 ./ mag_v2) * 180 / pi); % The first angle computed was for the second vertex, and the % last was for the first vertex. Scroll one position down to % make the last vertex be the first. angles = circshift(angles, [1, 0]); % Now determine if any vertices are concave and adjust the angles % accordingly. sgn = convex_angle_test(xy); % Any element of sgn that's -1 indicates that the angle is % concave. The corresponding angles have to be subtracted % from 360. I = find(sgn == -1); angles(I) = 360 - angles(I); %-------------------------------------------------------------------% function sgn = convex_angle_test(xy) % The rows of array xy are ordered vertices of a polygon. If the % kth angle is convex (>0 and <= 180 degress) then sgn(k) = % 1. Otherwise sgn(k) = -1. This function assumes that the first % vertex in the list is convex, and that no other vertex has a % smaller value of x-coordinate. These two conditions are true in % the first vertex generated by the MPP algorithm. Also the % vertices are assumed to be ordered in a clockwise sequence, and % there can be no duplicate vertices. % % The test is based on the fact that every convex vertex is on the % positive side of the line passing through the two vertices % immediately following each vertex being considered. If a vertex % is concave then it lies on the negative side of the line joining % the next two vertices. This property is true also if positive and % negative are interchanged in the preceding two sentences. % It is assumed that the polygon is closed. If not, close it. if size(xy, 1) == 1 | ~isequal(xy(1, :), xy(end, :)) xy(end + 1, :) = xy(1, :); end % Sign convention: sgn = 1 for convex vertices (i.e, interior angle > 0 % and <= 180 degrees), sgn = -1 for concave vertices. % Extreme points to be used in the following loop. A 1 is appended % to perform the inner (dot) product with w, which is 1-by-3 (see % below). L = 10^25; top_left = [-L, -L, 1]; top_right = [-L, L, 1]; bottom_left = [L, -L, 1]; bottom_right = [L, L, 1]; sgn = 1; % The first vertex is known to be convex. % Start following the vertices. for k = 2:length(xy) - 1 pfirst= xy(k - 1, :); psecond = xy(k, :); % This is the point tested for convexity. pthird = xy(k + 1, :); % Get the coefficients of the line (polygon edge) passing % through pfirst and psecond. w = polyedge(pfirst, psecond); % Establish the positive side of the line w1x + w2y + w3 = 0. % The positive side of the line should be in the right side of the % vector (psecond - pfirst). deltax and deltay of this vector % give the direction of travel. This establishes which of the % extreme points (see above) should be on the + side. If that % point is on the negative side of the line, then w is replaced by -w. deltax = psecond(:, 1) - pfirst(:, 1); deltay = psecond(:, 2) - pfirst(:, 2); if deltax == 0 & deltay == 0 error('Data into convexity test is 0 or duplicated.') end if deltax <= 0 & deltay >= 0 % Bottom_right should be on + side. vector_product = dot(w, bottom_right); % Inner product. w = sign(vector_product)*w; elseif deltax <= 0 & deltay <= 0 % Top_right should be on + side. vector_product = dot(w, top_right); w = sign(vector_product)*w; elseif deltax >= 0 & deltay <= 0 % Top_left should be on + side. vector_product = dot(w, top_left); w = sign(vector_product)*w; else % deltax >= 0 & deltay >= 0, so bottom_left should be on + side. vector_product = dot(w, bottom_left); w = sign(vector_product)*w; end % For the vertex at psecond to be convex, pthird has to be on the % positive side of the line. sgn(k) = 1; if (w(1)*pthird(:, 1) + w(2)*pthird(:, 2) + w(3)) < 0 sgn(k) = -1; end end %-------------------------------------------------------------------% function w = polyedge(p1, p2) % Outputs the coefficients of the line passing through p1 and % p2. The line is of the form w1x + w2y + w3 = 0. x1 = p1(:, 1); y1 = p1(:, 2); x2 = p2(:, 1); y2 = p2(:, 2); if x1==x2 w2 = 0; w1 = -1/x1; w3 = 1; elseif y1==y2 w1 = 0; w2 = -1/y1; w3 = 1; elseif x1 == y1 & x2 == y2 w1 = 1; w2 = 1; w3 = 0; else w1 = (y1 - y2)/(x1*(y2 - y1) - y1*(x2 - x1) + eps); w2 = -w1*(x2 - x1)/(y2 - y1); w3 = 1; end w = [w1, w2, w3]; %-------------------------------------------------------------------% function [xg, yg] = dupgone(x, y) % Eliminates duplicate, adjacent rows in [x y], except that the % first and last rows can be equal so that the polygon is closed. xg = x; yg = y; if size(xg, 1) > 2 I = find((x(1:end-1, :) == x(2:end, :)) & ... (y(1:end-1, :) == y(2:end, :))); xg(I) = []; yg(I) = []; end