Codebase list octave-iso2mesh / debian/1.9.6+ds-3 surfinterior.m
debian/1.9.6+ds-3

Tree @debian/1.9.6+ds-3 (Download .tar.gz)

surfinterior.m @debian/1.9.6+ds-3 — raw · history · blame 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function [pt,p0,v0,t,idx]=surfinterior(node,face)
%
% [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)
%

pt=[];
len=size(face,1);
for i=1:len
   p0=mean(node(face(i,1:3),:));
   plane=surfplane(node,face(i,:));
   v0=plane(1:3);

   [t,u,v]=raytrace(p0,v0,node,face(:,1:3));

   idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
   [ts, uidx]=unique(sort(t(idx)));
   if(~isempty(ts) && mod(length(ts),2)==0)
       ts=reshape(ts,[2 length(ts)/2]);
       tdiff=ts(2,:)-ts(1,:);
       [maxv,maxi]=max(tdiff);
       pt=p0+v0*(ts(1,maxi)+ts(2,maxi))*0.5;
       idx=idx(uidx);
       t=t(idx);
       break;
   end
end