function [R] = bqcurv(R, n),
%
% [R] = bqcurv(R, n) calculates the curvature of a mesh [R] at vertices
% [R.points] by fitting a bivariate quadratic polynomial to an n-ring
% around each vertex, where n is equal to one or two. If R contains the
% field "neighbours", a cell array the same length as R.points specifying
% the neighbours of each point
% R.points is a list of vertex co-ordinates in the form:
% pt. no., x-value, y-value, z-value, norm_x, norm_y, norm_z
% norm_* are components of the unit normal of that vertex
% and R.polygons is a list of vertices of the mesh polygons in the form:
% poly. no., pt. #1, pt. #2, pt. #3, Nx, Ny, Nz
% N is the unit normal of that polygon, norm being that of the vertex
%
% $Creation Date: 05/10/2004 $
% $Last Modified: 15/08/2006 $
% $Author: David Ferrari $
% $Contact: david_dot_ferrari_at_anu_dot_edu_dot_au
kappa_max = zeros(size(R.points,1),1);
kappa_min = zeros(size(R.points,1),1);
wupdate = round(length(R.points)/100); % Update counter for the waitbar
h = waitbar(0,'Computing Curvatures...');
for pt1 = 1:size(R.points,1), %for each vertex
if mod(pt1,wupdate)==0,
% Update waitbar every wupdate'th time
waitbar(pt1/length(R.points));
end %if mod(poly,wupdate)
% Establish the coordinates of this vertex, and its normal
x_0 = R.points(pt1, 2);
y_0 = R.points(pt1, 3);
z_0 = R.points(pt1, 4);
Nx = R.points(pt1,5);
Ny = R.points(pt1,6);
Nz = R.points(pt1,7);
%#######################################################################
% 1) Create a list of neighbouring points (1-ring or 2-ring)
if isfield(R, 'neighbours'),
%if a list of neighbours has already been created
neighbours = R.neighbours{pt1};
else %create a list of neighbours
neighbours = find_neighbours(pt1, R.polygons, n);
R.neighbours{pt1} = neighbours;
end %if isfield(R, 'neighbours'),
if (size(neighbours,1)) <6, % if pt1 is of valency < 6
neighbours = find_neighbours(neighbours, R.polygons, 2);
end
%#######################################################################
% 2) Compute the (vertex) normal N (unnecessary - already defined in
% .wrl format)
%#######################################################################
% 3) Compute coefficients of the "pseudo-tangent" plane P (normal to N)
% with the plane defined implicitly as:
% n.(x-y_0) = Nx(x-x_0) + Ny(y-y_0) + Nz(z-z_0)
% = Nx.x + Ny.y + Nz.z - (Nx.x_0 + Ny.y_0 + Nz.z_0)
% = Ax + By + Cz + D = 0
A = Nx; B = Ny; C = Nz;
D = -(Nx*x_0 + Ny*y_0 + Nz*z_0);
%#######################################################################
% 4) Define an arbitrary orthonormal coordinate system in P
% with pt1 as its origin
if Nx ~= 0, %if the x-component of the normal is non-zero
a = 1/Nx*[-(Ny+Nz), Nx, Nx];
elseif Ny ~= 0,
a = 1/Ny*[Ny, -(Nx+Nz), Ny];
elseif Nz ~= 0,
a = 1/Nz*[Nz, Nz, -(Nx+Ny)];
end %if Nx ~= 0, %if the x-component of the normal is non-zero
% Compute the basis vectors on the pseudo-tangent plane
b1 = a/norm(a);
b2 = cross([Nx,Ny,Nz],b1);
%#######################################################################
% 5) Compute the distances of all vertices in the neighbourhood from P
dist = zeros(size(neighbours,1),1); % Initilise the distance array
for i = 1:size(neighbours,1), % For each point in the neighbourhood
x_i = R.points(neighbours(i), 2);
y_i = R.points(neighbours(i), 3);
z_i = R.points(neighbours(i), 4);
% Calculate the distance
dist(i) = A*x_i + B*y_i + C*z_i + D;
end %for i = 1:size(neighbours,1),
%#######################################################################
% 6) Project all points in the neighbourhood onto P, and represent the
% coordinates of this projection in the local coordinate system defined
% in step (3)
pp = zeros(size(neighbours,1),3); % Initialise the projected point array
for i = 1:size(neighbours,1), % For each point in the neighbourhood
x_i = R.points(neighbours(i), 2);
y_i = R.points(neighbours(i), 3);
z_i = R.points(neighbours(i), 4);
pp(i,1) = x_i - dist(i)*Nx; % The x-value of the projected point
pp(i,2) = y_i - dist(i)*Ny; % The y-value of the projected point
pp(i,3) = z_i - dist(i)*Nz; % The z-value of the projected point
end %for i = 1:size(neighbours,1),
%#######################################################################
% 7) Interpret projections of each point in P as abscissa values, and
% the distance of each point from P as ordinate values (di)
% We already have the ordinates in (dist), so we calculate the abscissa.
abscissa = zeros(size(neighbours,1),2); % Initialise the abscissa array
for i = 1:size(neighbours,1), % For each point in the neighbourhood
% The difference vector between this projected point and the origin
diff = pp(i,:) - [x_0,y_0,z_0];
% The abscissae as defined by the new coordinate system on P:
abscissa(i,1) = dot(diff, b1);
abscissa(i,2) = dot(diff, b2);
end %for i = 1:size(neighbours,1),
%#######################################################################
% 8) Construct a bivariate quadratic polynomial by solving the least
% squares problem Uc = dist
U = zeros(size(neighbours,1),3); % Initialise the "U" matrix
for i = 1:size(neighbours,1), % For each point in the neighbourhood
u_i = abscissa(i,1);
v_i = abscissa(i,2);
U(i,:) = [u_i^2, u_i*v_i, v_i^2];
end %for i = 1:size(U,1),
% We now solve:
% U*c = dist, ie
% transpose(U)*U*c = transpose(U)*(dist)
c = inv(transpose(U)*U) * transpose(U)*dist;
%XXXXXXXXXXXXXXXXXX transpose(U)*U might be simgular XXXXXXXXXXXXXXXXXXX
% We now have the polynomial f:
% f(u,v) = 1/2*(c(1)*u^2 + 2*c(2)*u*v + c(3)*v^2)
%#######################################################################
% 9) Compute the curvatures from the bivariate polynomial
% The principle curvatures are given the two real roots of
% k^2 - (c(1) + c(3))*k - c(2)^2 = 0
if any(isnan(c)) || any(isinf(c)), % If U was singular
kappa_max(pt1) = 0;
kappa_min(pt1) = 0;
else
kappa_max(pt1) = ((c(1)+c(3)) + sqrt((c(1)+c(3))^2-4*1*(c(1)*c(3)...
-c(2)^2)))/2;
kappa_min(pt1) = ((c(1)+c(3)) - sqrt((c(1)+c(3))^2-4*1*(c(1)*c(3)...
-c(2)^2)))/2;
end %if any(isnan(c)) | any(isinf(c)),
end %for pt1 = 1:size(R.points,1), %for each vertex
close(h);
R.kappa_max = kappa_max;
R.kappa_min = kappa_min;
end
%#######################################################################
function neighbours = find_neighbours(pts, polygons, n),
% neighbours = find_neighbours(pts, polygons, n) creates a list of "mesh"
% neighbours in the n-ring of the point list pts, where polygons is of the
% form:
% poly. no., pt. #1, pt. #2, pt. #3, Nx, Ny, Nz
% N is the unit normal of that polygon, norm being that of the vertex
%
% $Creation Date: 5/10/2004 $
% $Last Modified: 14/10/2004 $
% $Author: David Ferrari $
%#######################################################################
% 1) Create a list of neighbouring points (1-ring or 2-ring)
if n==1, % if a 1-ring neighbourhood required
A = ismember(polygons(:,2:4), pts);
[polynums, temp] = find(A);
% polynums is a list of polygon no's with pts as vertices
B = ~ismember(polygons(polynums,2:4),pts);
neighbours = B.*polygons(polynums,2:4);
elseif n>1, % if a larger-ring neighbourhood required
% Find all neighbours 'pts2' in a 1-ring around all pts
A = ismember(polygons(:,2:4), pts);
[polynums, temp] = find(A);
% polynums is a list of polygon no's with pts as vertices
B = ~ismember(polygons(polynums,2:4),pts);
pts2 = B.*polygons(polynums,2:4);
pts2 = unique(pts2);
pts2 = pts2(:);
% Recursively call the function find_neighbours
neighbours = find_neighbours(pts2, polygons, n-1);
neighbours = [neighbours;pts2];
end %if n==1, ...
neighbours = unique(neighbours(:));
if neighbours(1) == 0,
neighbours = neighbours(2:size(neighbours,1)); %remove zeros
end
% Remove the original list pts from the list neighbours
already_counted = ismember(neighbours, pts);
neighbours(already_counted) = [];