-
Notifications
You must be signed in to change notification settings - Fork 326
/
compute_normal.sci
57 lines (48 loc) · 1.65 KB
/
compute_normal.sci
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
function [normal,normalf] = compute_normal(vertex,face)
// compute_normal - compute the normal of a triangulation
//
// [normal,normalf] = compute_normal(vertex,face);
//
// normal(i,:) is the normal at vertex i.
// normalf(j,:) is the normal at face j.
//
// Copyright (c) 2004 Gabriel Peyré
[vertex,face] = check_face_vertex(vertex,face);
nface = size(face,2);
nvert = size(vertex,2);
normal = zeros(3, nvert);
// unit normals to the faces
normalf = crossp( vertex(:,face(2,:))-vertex(:,face(1,:)), ...
vertex(:,face(3,:))-vertex(:,face(1,:)) );
d = sqrt( sum(normalf.^2,1) ); d(d<1e-10)=1;
normalf = normalf ./ mtlb_repmat( d, 3,1 );
// unit normal to the vertex
normal = zeros(3,nvert);
for i=1:nface
f = face(:,i);
for j=1:3
normal(:,f(j)) = normal(:,f(j)) + normalf(:,i);
end
end
// normalize
d = sqrt( sum(normal.^2,1) ); d(d<1e-10)=1;
normal = normal ./ mtlb_repmat( d, 3,1 );
// enforce that the normal are outward
v = vertex - mtlb_repmat(mean(vertex,1), 3,1);
s = sum( v.*normal, 2 );
if sum(s>0)<sum(s<0)
// flip
normal = -normal;
normalf = -normalf;
end
endfunction
//////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////
function z = crossp(x,y)
// x and y are (m,3) dimensional
z = x;
z(1,:) = x(2,:).*y(3,:) - x(3,:).*y(2,:);
z(2,:) = x(3,:).*y(1,:) - x(1,:).*y(3,:);
z(3,:) = x(1,:).*y(2,:) - x(2,:).*y(1,:);
endfunction