-
Notifications
You must be signed in to change notification settings - Fork 326
/
check_incircle.m
50 lines (46 loc) · 1.39 KB
/
check_incircle.m
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
function ic = check_incircle(vertex,face,points)
% check_incicle_edge - compute "empty circle" property for a set of faces
%
% ic = check_incicle(vertex,points);
% -> ic(i)==1 if the circumcircle of the 3 points
% vertex(:,face(:,i)) does not contain the point vertex(:,points(i))
%
% Copyright (c) 2008 Gabriel Peyre
n = size(face, 2);
if n~=length(points)
error('Problem');
end
if 0
%% Apparently there is a bug here in the code %%
% set of points
a1 = vertex(:,face(1,:));
a2 = vertex(:,face(2,:));
a3 = vertex(:,face(3,:));
p = vertex(:,points);
% build matrix for determinant
if 1
A = zeros(3,3,n);
A(1:2,1,:) = reshape(a1-p, [2 1 n]);
A(1:2,2,:) = reshape(a2-p, [2 1 n]);
A(1:2,3,:) = reshape(a3-p, [2 1 n]);
d = [sum((a1-p).^2); sum((a2-p).^2); sum((a3-p).^2)];
A(3,:,:) = reshape( d, [1 3 n] );
ic = det3(A)>0;
else
A = ones(4,4,n);
A(1:2,1,:) = reshape(a1, [2 1 n]);
A(1:2,2,:) = reshape(a2, [2 1 n]);
A(1:2,3,:) = reshape(a3, [2 1 n]);
A(1:2,4,:) = reshape(p, [2 1 n]);
d = [sum(a1.^2); sum(a2.^2); sum(a3.^2); sum(p.^2)];
A(4,:,:) = reshape( d, [1 4 n] );
ic = det3(A)<0;
end
end
% orthocenter
[q,r] = compute_orthocenter(vertex,face);
% distance to centerpoint
P = vertex(:,points);
dc = sum( (q-P).^2, 1);
% check in/out
ic = dc>r.^2;