surfinterior

 [pt,p0,v0,t,idx]=surfinterior(node,face)

 identify a point that is enclosed by the (closed) surface

 author: Qianqian Fang, <q.fang at neu.edu>

 input:
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)

 output:
   pt: the interior point coordinates [x y z]
   p0: ray origin used to determine the interior point
   v0: the vector used to determine the interior point
   t : ray-tracing intersection distances (with signs) from p0. the
       intersection coordinates can be expressed as p0+t(i)*v0
   idx: index to the face elements that intersect with the ray, order
       match that of t

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

0001 function [pt,p0,v0,t,idx]=surfinterior(node,face)
0002 %
0003 % [pt,p0,v0,t,idx]=surfinterior(node,face)
0004 %
0005 % identify a point that is enclosed by the (closed) surface
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 %
0009 % input:
0010 %   node: a list of node coordinates (nn x 3)
0011 %   face: a surface mesh triangle list (ne x 3)
0012 %
0013 % output:
0014 %   pt: the interior point coordinates [x y z]
0015 %   p0: ray origin used to determine the interior point
0016 %   v0: the vector used to determine the interior point
0017 %   t : ray-tracing intersection distances (with signs) from p0. the
0018 %       intersection coordinates can be expressed as p0+t(i)*v0
0019 %   idx: index to the face elements that intersect with the ray, order
0020 %       match that of t
0021 %
0022 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0023 %
0024 
0025 pt=[];
0026 len=size(face,1);
0027 for i=1:len
0028    p0=mean(node(face(i,1:3),:));
0029    plane=surfplane(node,face(i,:));
0030    v0=plane(1:3);
0031 
0032    [t,u,v]=raytrace(p0,v0,node,face(:,1:3));
0033 
0034    idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
0035    [ts, uidx]=unique(sort(t(idx)));
0036    if(~isempty(ts) && mod(length(ts),2)==0)
0037        ts=reshape(ts,[2 length(ts)/2]);
0038        tdiff=ts(2,:)-ts(1,:);
0039        [maxv,maxi]=max(tdiff);
0040        pt=p0+v0*(ts(1,maxi)+ts(2,maxi))*0.5;
0041        idx=idx(uidx);
0042        t=t(idx);
0043        break;
0044    end
0045 end