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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
function [newpt elemid weight] = proj2mesh(v, f, pt, nv, cn, radmax)
% [newpt elemid weight]=proj2mesh(v,f,pt,nv,cn)
%
% project a point cloud on to the surface mesh (surface can only be triangular)
%
% author: Qianqian Fang <q.fang at neu.edu>
% date: 12/12/2008
%
% parameters:
% v: node coordinate of the surface mesh (nn x 3)
% f: element list of the surface mesh (3 columns for
% triangular mesh, 4 columns for cubic surface mesh)
% pt: points to be projected, 3 columns for x,y and z respectively
% nv: nodal norms (vector) calculated from nodesurfnorm.m
% with dimensions of (size(v,1),3)
% cn: a integer vector with the length of p, denoting the closest
% surface nodes (indices of v) for each point in p. this
% value can be calculated from dist2surf.m
% radmax: if speicified, the search for elements to project will be
% limited to those within a bounding box with half-edge-length
% of radmax centered at the point to be projected
%
% if nv and cn are not supplied, proj2mesh will project the point
% cloud onto the surface by the direction pointing to the centroid
% of the mesh
%
% outputs:
% newpt: the projected points from p
% elemid: a vector of length of p, denotes which surface trangle (in elem)
% contains the projected point
% weight: the barycentric coordinate for each projected points, these are
% the weights
%
% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
% this function is part of "metch" toobox, see COPYING for license
cent = mean(v);
enum = length(f);
ec = reshape(v(f(:, 1:3)', :)', [3 3, enum]);
centroid = squeeze(mean(ec, 2));
newpt = zeros(size(pt, 1), 3);
elemid = zeros(size(pt, 1), 1);
weight = zeros(size(pt, 1), 3);
[idoldmesh, loc] = ismember(pt, v, 'rows');
idnode = find(idoldmesh);
if (~isempty(idnode))
[tt, ll] = ismember(loc(idnode), f);
[p1, p2] = ind2sub(size(f), ll); % p1 is the index in f
newpt(idnode, :) = pt(idnode, :);
elemid(idnode) = p1;
weight(sub2ind(size(weight), idnode, p2)) = 1;
end
radlimit = -1;
if (nargin >= 5)
% if nv and cn are supplied, use nodal norms to project the points
direction = nv(cn, :);
if (nargin >= 6)
radlimit = radmax;
end
elseif (nargin == 3)
% otherwise, project toward the centroid
direction = pt - repmat(cent, size(pt, 1), 1);
end
for t = 1:size(pt, 1)
if (idoldmesh(t) ~= 0)
continue
end
maxdist = sqrt(sum((pt(t, :) - cent) .* (pt(t, :) - cent)));
if (radlimit > 0)
maxdist = radlimit;
end
idx = find(centroid(1, :) > pt(t, 1) - maxdist & ...
centroid(1, :) < pt(t, 1) + maxdist & ...
centroid(2, :) > pt(t, 2) - maxdist & ...
centroid(2, :) < pt(t, 2) + maxdist & ...
centroid(3, :) > pt(t, 3) - maxdist & ...
centroid(3, :) < pt(t, 3) + maxdist);
dist = centroid(:, idx);
dist(1, :) = dist(1, :) - pt(t, 1);
dist(2, :) = dist(2, :) - pt(t, 2);
dist(3, :) = dist(3, :) - pt(t, 3);
c0 = sum(dist .* dist);
% sort the distances to accelate the calculation
[c1, sorted] = sort(c0);
for i = 1:length(idx)
% project the point along vector direction and calculate the intersection to a plane
[inside, p, w] = linextriangle(pt(t, :), pt(t, :) + direction(t, :), v(f(idx(sorted(i)), :), :));
% the intersection is within the current trianglar surface
if (inside)
newpt(t, :) = p;
weight(t, :) = w;
elemid(t) = idx(sorted(i));
break
end
end
end
|
