internalpoint

 p=internalpoint(v,aloop)

 imperical function to find an internal point
 of a planar polygon

 author: Qianqian Fang, <q.fang at neu.edu>
 date: 2008/04/08

 input:   
    v:     x,y,z coordinates of each node of the mesh
    aloop:  input, a single vector separated by NaN, each segment
             is a close-polygon consisted by node IDs 
 output:
    p:   output, [x y z] of an internal point of aloop

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function p=internalpoint(v,aloop)
0002 %
0003 % p=internalpoint(v,aloop)
0004 %
0005 % imperical function to find an internal point
0006 % of a planar polygon
0007 %
0008 % author: Qianqian Fang, <q.fang at neu.edu>
0009 % date: 2008/04/08
0010 %
0011 % input:
0012 %    v:     x,y,z coordinates of each node of the mesh
0013 %    aloop:  input, a single vector separated by NaN, each segment
0014 %             is a close-polygon consisted by node IDs
0015 % output:
0016 %    p:   output, [x y z] of an internal point of aloop
0017 %
0018 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0019 %
0020 
0021 p=[];
0022 nd=v(aloop,:);
0023 boxfacet=find(sum(abs(diff(nd)))<1e-2); % find a flat loop
0024 if(length(boxfacet))   % if the loop is flat along x/y/z dir
0025      bf=boxfacet(1);    % no degeneracy allowed
0026      idx=setdiff([1 2 3],bf);
0027      
0028      p0=(nd(1,:)+nd(2,:))/2;
0029      pvec=complex(p0(idx(1)),p0(idx(2)));
0030      vec=nd(2,:)-nd(1,:);
0031      vec=complex(vec(idx(1)),vec(idx(2)))*exp(i*pi/2)*(1e-5)/sqrt(sum(vec.*vec));
0032      testpt=[real(pvec+vec) imag(pvec+vec);real(pvec-vec) imag(pvec-vec)];
0033      in=inpolygon(testpt(:,1),testpt(:,2), nd(:,idx(1)),nd(:,idx(2)));
0034      p=testpt(find(in>0),:);
0035      p([bf,idx(1),idx(2)])=[nd(1,bf),p];
0036 end
0037 
0038 if(length(p)==0|length(p)==2) 
0039     error('fail to find an internal point of curve');
0040 end