matlab_class/gcmfaces_misc/gcmfaces_polygonangle.m

function [angsum]=gcmfaces_polygonangle(px,py,x,y);
%[angsum]=gcmfaces_polygonangle(px,py,x,y);
%object:    compute sum of interior angles for polygons (when input
%           is px,py) or points vs polygons (when input is px,py,x,y)
%inputs:    px,py are MxN matrices where each line specifies one polygon
%(optional) x,y are position vectors
%outputs:   ang are the corresponding sums of interior angles

M=size(px,1); N=size(px,2);
doPointsInPolygon=0; P=1;
if nargin>2;
    doPointsInPolygon=1;
    sizxy=size(x);
    x=reshape(x,[1 prod(sizxy)]);
    y=reshape(y,[1 prod(sizxy)]);
    P=length(x);
end;

angsum=zeros(M,P);
for ii=0:N-1;
    ppx=circshift(px,[0 -ii]);
    ppy=circshift(py,[0 -ii]);
    
    if ~doPointsInPolygon;
        %compute sum of corner angles
    v1x=ppx(:,2)-ppx(:,1); 
    v1y=ppy(:,2)-ppy(:,1);
    v2x=ppx(:,4)-ppx(:,1); 
    v2y=ppy(:,4)-ppy(:,1);
    else;
        %compute sum of sector angles        
        v1x=ppx(:,1)*ones(1,P)-ones(M,1)*x; 
        v1y=ppy(:,1)*ones(1,P)-ones(M,1)*y;
        v2x=ppx(:,2)*ones(1,P)-ones(M,1)*x; 
        v2y=ppy(:,2)*ones(1,P)-ones(M,1)*y;
    end;
    g_acos=acos( ( v1x.*v2x+v1y.*v2y )./sqrt( v1x.*v1x+v1y.*v1y )./sqrt( v2x.*v2x+v2y.*v2y ) );
    g_sin= ( v1x.*v2y-v1y.*v2x )./sqrt( v1x.*v1x+v1y.*v1y )./sqrt( v2x.*v2x+v2y.*v2y ) ;
    angsum=angsum+rad2deg(g_acos.*sign(g_sin));    
end;

1	gforget	1.1	function [angsum]=gcmfaces_polygonangle(px,py,x,y);
2			%[angsum]=gcmfaces_polygonangle(px,py,x,y);
3			%object: compute sum of interior angles for polygons (when input
4			% is px,py) or points vs polygons (when input is px,py,x,y)
5			%inputs: px,py are MxN matrices where each line specifies one polygon
6			%(optional) x,y are position vectors
7			%outputs: ang are the corresponding sums of interior angles
8
9			M=size(px,1); N=size(px,2);
10			doPointsInPolygon=0; P=1;
11			if nargin>2;
12			doPointsInPolygon=1;
13			sizxy=size(x);
14			x=reshape(x,[1 prod(sizxy)]);
15			y=reshape(y,[1 prod(sizxy)]);
16			P=length(x);
17			end;
18
19			angsum=zeros(M,P);
20			for ii=0:N-1;
21			ppx=circshift(px,[0 -ii]);
22			ppy=circshift(py,[0 -ii]);
23
24			if ~doPointsInPolygon;
25			%compute sum of corner angles
26			v1x=ppx(:,2)-ppx(:,1);
27			v1y=ppy(:,2)-ppy(:,1);
28			v2x=ppx(:,4)-ppx(:,1);
29			v2y=ppy(:,4)-ppy(:,1);
30			else;
31			%compute sum of sector angles
32			v1x=ppx(:,1)ones(1,P)-ones(M,1)x;
33			v1y=ppy(:,1)ones(1,P)-ones(M,1)y;
34			v2x=ppx(:,2)ones(1,P)-ones(M,1)x;
35			v2y=ppy(:,2)ones(1,P)-ones(M,1)y;
36			end;
37			g_acos=acos( ( v1x.v2x+v1y.v2y )./sqrt( v1x.v1x+v1y.v1y )./sqrt( v2x.v2x+v2y.v2y ) );
38			g_sin= ( v1x.v2y-v1y.v2x )./sqrt( v1x.v1x+v1y.v1y )./sqrt( v2x.v2x+v2y.v2y ) ;
39			angsum=angsum+rad2deg(g_acos.*sign(g_sin));
40			end;
41