-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtricurv_v01.m
199 lines (177 loc) · 6.09 KB
/
tricurv_v01.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
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
function out=tricurv_v01(tri,p)
% Function to calculate the principal curvatures, and their respective
% directions on a triangular mesh. Approximations of curvature are based on
% local (N=1) neighborhood elements and vertices.
% Note that calculations at vertices with few adjacent triangles, and hence
% few adjacent vertices, are expanded to a greater neighborhood.
%
%
% Reference:
% 1) Chen and Schmitt (1992) Intrinsic surface properties from surface
% triangulation
% 2) Dong et al. (2005) Curvature estimation on triangular mesh,
% JZUS
%
% This code makes use of routines: buildInverseTriangualtion.m & removeDO.m
% initially written by: David Gringas. He is gratefully acknowledged
%
% Input: t <mx3> array of triangle indices into:
% xyz <nx3> coordinates
% Output: structure containing: Principal curvatures and directions,
% Chen and Schmitt's coefficients, etc.
% Version: 1
% JOK 030809
% I/O check:
if nargin~=2
error('Wrong # of input')
end
if nargout ~= 1
error('Output is a strcture, wrong designation!')
end
nt= size(tri);
if nt(2)~=3
error('Triangle element matrix should be mx3!')
end
mp = size(p);
if mp(2)~=3
error('Vertices should be nx3!')
end
% 1) Compute vertex normals (weighted by distance)
nv=vertnorm_dw(tri,p);
% 2) Define neighborhoods and calculate weighted normals
% Average at vertices
itri=buildInverseTriangulation(tri);
for j=1:length(p)
ind01=removeD0(itri(j,:));% Triangles adjacent to vertex j
ind02 = tri(ind01,:);ind02 = unique(ind02(:));% Vertices of adjacent triangles
if length(ind02)<5 % Larger neighborhood for triangles with less than 5 unique triangles
for i=1:length(ind02)
indaux01 = removeD0(itri(ind02(i),:));% Triangles adjacent to neighborhood vertices in ind02(i)
indaux02 = tri(indaux01,:);% Vertices of larger neighborhood
indaux02 = unique(indaux02(:));
ind02 = unique([ind02;indaux02]);
end
end
% Define tangent vectors
pdist = p(ind02,:)-repmat(p(j,:),length(ind02),1);
aux01 = pdist-repmat(dot(pdist,nv(ind02,:),2),1,3).*nv(ind02,:);
t = aux01./repmat(sqrt(sum(aux01.*aux01,2)),1,3);
% Normal curvature
kn = -dot(pdist,nv(ind02,:)- ...
repmat(nv(j,:),length(ind02),1),2)./dot(pdist,pdist,2);
[mkn,ind03] = max(kn);
% Local principal directions
out.e1(j,:) = t(ind03,:);
aux02 = cross(out.e1(j,:),nv(j,:),2);
out.e2(j,:) = aux02./repmat(sqrt(sum(aux02.*aux02,2)),1,3);
% Chen and Schmitt: calculation of coefficients
thetai = real(acos(dot(t,repmat(out.e1(j,:),length(ind02),1),2)));
% Check for NaN
thetai(isnan(thetai)) = 0; % Can be set to zero as all coefficients
% contain a sin term and thus are zero
kn(isnan(kn))=0;
a11(j,1) = sum((cos(thetai)).^2.*(sin(thetai)).^2,1);
a12(j,1) = sum(cos(thetai).*(sin(thetai)).^3,1);
a22(j,1) = sum((sin(thetai)).^4,1);
a13(j,1) = sum((kn-kn(ind03)*(cos(thetai)).^2).*cos(thetai).*sin(thetai),1);
a23(j,1) = sum((kn-kn(ind03)*(cos(thetai)).^2).*(sin(thetai)).^2,1);
% coefficients
a(j,1) = kn(ind03);
b(j,1) = (a13(j)*a22(j)-a23(j)*a12(j))/(a11(j)*a22(j)-a12(j)^2);
c(j,1) = (a11(j)*a23(j)-a12(j)*a13(j))/(a11(j)*a22(j)-a12(j)^2);
end
% Gaussian, mean normal, principal curvatures
aux03 = sort([.5*(a+c+sqrt((a-c).^2+4*b.^2)),.5*(a+c-sqrt((a-c).^2+4*b.^2))],2);
out.k1 = aux03(:,1);
out.k2 = aux03(:,2);
out.km = .5*(out.k1+out.k2);
out.kg = out.k1.*out.k2;
end % End tricurv_v01
%%%%%%%%%%%%%% SUBFUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunctions
function invTRI=buildInverseTriangulation(TRI)
% Building the inverse triangulation, i.e. a link from node indexes to
% triangle indexes.
nbTri=length(TRI);
nbNode=max(reshape(TRI,[],1));
comp=zeros(nbNode,1);
invTRI=zeros(nbNode,8);
for i=1:nbTri
for j=1:3
index=TRI(i,j);
comp(index)=comp(index)+1;
invTRI(index,comp(index))=i;
end
end
end % End buildInverseTriangulation
%%
function out=removeD0(x)
% Removing duplicate and null values
s=sort(x);
s1=[0,s];
s2=[s,s(length(s))];
ds=(s1-s2);
out=s(logical(ds~=0));
end% End removeD0
%%
function nv=vertnorm_dw(tri,p)
% Function to compute normal vector of vertices comprising a triangular
% mesh. Based on trinormal and computeNormalVectorTriangulation.m by David
% Gringas
% Input: tri mx3 <triangular index matrix>
% p nx3 array of vertices
% Output: nvec <nx3> array of normal vectors
% JOK230709
% Version: 1
% I/O check
% to be completed
% Construct vectors
v = [p(tri(:,3),1)-p(tri(:,1),1), p(tri(:,3),2)-p(tri(:,1),2), p(tri(:,3),3)-p(tri(:,1),3)];
w = [p(tri(:,2),1)-p(tri(:,1),1), p(tri(:,2),2)-p(tri(:,1),2), p(tri(:,2),3)-p(tri(:,1),3)];
% Calculate cross product
normvec = [v(:,2).*w(:,3)-v(:,3).*w(:,2), ...
-(v(:,1).*w(:,3)-v(:,3).*w(:,1)), ...
v(:,1).*w(:,2)-v(:,2).*w(:,1)];
% Normalize
lnvec = sqrt(sum(normvec.*normvec,2));
nvec = normvec./repmat(lnvec,1,3);
% Average at vertices
itri=buildInverseTriangulation(tri);
nvecx=zeros(length(p),1);
nvecy=zeros(length(p),1);
nvecz=zeros(length(p),1);
for j=1:length(p)
% Find centroids and weight based on distance:
ind01=removeD0(itri(j,:));
cen=tricentroid(p,tri(ind01,:));
nc=size(cen);
distcen = cen-repmat(p(j,:),nc(1),1);
w = 1./sqrt(sum(distcen.*distcen,2));
nvecx=mean(w.*nvec(ind01,1));
nvecy=mean(w.*nvec(ind01,2));
nvecz=mean(w.*nvec(ind01,3));
% Nomalize
lnv = sqrt(nvecx^2+nvecy^2+nvecz^2);
nv(j,:) = [nvecx/lnv,nvecy/lnv,nvecz/lnv];
end
end % End vertnorm
%%
function out = tricentroid(v,tri)
% Function to output the centroid of triangluar elements.
% Note that the output will be of length(tri)x3
% Input: <v> nx2 or 3: vertices referenced in tri
% <tri> mx3: triangle indices
% Version: 1
% JOK 300509
% I/O check
[nv,mv]=size(v);
[nt,mt]=size(tri);
if mv==2
v(:,3) = zeros(nv,1);
elseif mt~=3
tri=tri';
end
out(:,1) = 1/3*(v(tri(:,1),1)+v(tri(:,2),1)+v(tri(:,3),1));
out(:,2) = 1/3*(v(tri(:,1),2)+v(tri(:,2),2)+v(tri(:,3),2));
out(:,3) = 1/3*(v(tri(:,1),3)+v(tri(:,2),3)+v(tri(:,3),3));
end%tricentroid