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
|
function p = internalpoint(v, aloop)
%
% 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)
%
p = [];
nd = v(aloop, :);
boxfacet = find(sum(abs(diff(nd))) < 1e-2); % find a flat loop
if (~isempty(boxfacet)) % if the loop is flat along x/y/z dir
bf = boxfacet(1); % no degeneracy allowed
idx = setdiff([1 2 3], bf);
p0 = (nd(1, :) + nd(2, :)) / 2;
pvec = complex(p0(idx(1)), p0(idx(2)));
vec = nd(2, :) - nd(1, :);
vec = complex(vec(idx(1)), vec(idx(2))) * exp(1i * pi / 2) * (1e-5) / sqrt(sum(vec .* vec));
testpt = [real(pvec + vec) imag(pvec + vec); real(pvec - vec) imag(pvec - vec)];
in = inpolygon(testpt(:, 1), testpt(:, 2), nd(:, idx(1)), nd(:, idx(2)));
p = testpt(in > 0, :);
p([bf, idx(1), idx(2)]) = [nd(1, bf), p];
end
if (isempty(p) || length(p) == 2)
error('fail to find an internal point of curve');
end
|
