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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
function [t, u, v, idx, xnode] = raysurf(p0, v0, node, face)
%
% [t,u,v,idx,xnode]=raysurf(p,v,node,face)
%
% perform a Havel-styled ray tracing for a triangular surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
% p0: list of starting points of the rays
% v0: directional vector of the rays,
% node: a list of node coordinates (nn x 3)
% face: a surface mesh triangle list (ne x 3)
%
% output:
% t: distance from p0 to the intersection point for each surface
% triangle, if t(i)=NaN, no intersection was found for that ray
% u: bary-centric coordinate 1 of all intersection points
% v: bary-centric coordinate 2 of all intersection points
% the final bary-centric triplet is [u,v,1-u-v]
% idx: idx lists the IDs of the face elements that intersects
% each ray
% xnode: optional output, if requested, xnode gives the intersection
% point coordinates; to compute manually, xnode=p0+repmat(t,1,3).*v0
%
% Reference:
% [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using
% SSE4)," IEEE Trans. on Visualization and Computer Graphics,
% 16(3):434-438 (2010)
% [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous
% Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
% press)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
len = size(p0, 1);
if (len == 0)
error('p0 can not be empty');
end
if (size(node, 2) < 3)
error('node must contain at least 3 columns');
end
if (size(face, 2) < 3)
error('face must contain at least 3 columns');
end
if (size(v0, 1) == 1 || size(v0, 2) == 1 && len > 1)
v0 = repmat(v0(:)', len, 1);
end
t = zeros(len, 1) * nan;
u = t;
v = t;
idx = t;
for i = 1:len
[ti, ui, vi, id] = raytrace(p0(i, :), v0(i, :), node, face);
if (isempty(id))
continue
end
ti = ti(id);
tpid = find(ti >= 0);
if (isempty(tpid))
continue
end
[tmin, tloc] = min(ti(find(ti >= 0)));
t(i) = tmin;
u(i) = ui(id(tpid(tloc)));
v(i) = vi(id(tpid(tloc)));
idx(i) = id(tpid(tloc));
end
if (nargout >= 5)
xnode = p0 + repmat(t, 1, 3) .* v0;
end
|
