spherical_conformal_map

Spherical Conformal Map of a closed genus-0 surface. Methodology details please refer to [1]. Some code is derived from Gary Choi.

  1. K.C. Lam, P.T. Choi and L.M. Lui, FLASH: Fast landmark aligned spherical harmonic parameterization for genus-0 closed brain surfaces, UCLA CAM Report, ftp://ftp.math.ucla.edu/pub/camreport/cam13-79.pdf?

Contents

Syntax

uvw = spherical_conformal_map(face,vertex)

Description

face  : double array, nf x 3, connectivity of mesh
vertex: double array, nv x 3, vertex of mesh
uvw: double array, nv x 3, spherical uvw coordinates of vertex on 3D unit sphere

Contribution

Author : Wen Cheng Feng
Created: 2014/03/27
Revised: 2014/03/28 by Wen, add doc
Copyright 2014 Computational Geometry Group
Department of Mathematics, CUHK
http://www.math.cuhk.edu.hk/~lmlui
function uvw = spherical_conformal_map(face,vertex)
dv1 = vertex(face(:,2),:) - vertex(face(:,3),:);
e1 = sqrt(dot(dv1,dv1,2));
dv2 = vertex(face(:,3),:) - vertex(face(:,1),:);
e2 = sqrt(dot(dv2,dv2,2));
dv3 = vertex(face(:,1),:) - vertex(face(:,2),:);
e3 = sqrt(dot(dv3,dv3,2));
e123 = e1+e2+e3;
regularity = abs(e1./e123-1/3)+abs(e2./e123-1/3)+abs(e3./e123-1/3);
% choose vertex bi as the big triangle
[~,bi] = min(regularity);
nv = size(vertex,1);
A = laplace_beltrami(face,vertex);
fi = face(bi,:);
[I,J,V] = find(A(fi,:));
A = A - sparse(fi(I),J,V,nv,nv) + sparse(fi,fi,[1,1,1],nv,nv);

% Set boundary condition for big triangle
x1 = 0; y1 = 0;
x2 = 100; y2 = 0;
a = vertex(fi(2),:) - vertex(fi(1),:);
b = vertex(fi(3),:) - vertex(fi(1),:);
ratio = norm([x1-x2,y1-y2])/norm(a);
y3 = norm([x1-x2,y1-y2])*norm(cross(a,b))/norm(a)^2;
x3 = sqrt(norm(b)^2*ratio^2-y3^2);

% Solve matrix equation
% d = complex(zeros(nv,1));
% d(fi) = [x1+1i*y1;x2+1i*y2;x3+1i*y3];
% z = A\d;
d = zeros(nv,2);
d(fi,:) = [x1,y1;x2,y2;x3,y3];
uv = A\d;
z = uv(:,1) + 1i*uv(:,2);
z = z - mean(z);
dz2 = abs(z).^2;
vertex_new = [2*real(z)./(1+dz2), 2*imag(z)./(1+dz2), (-1+dz2)./(1+dz2)];

% Find optimal big triangle size
% Reason: the distribution will be the best
% if the southmost triangle has similar size of the northmost one
w = complex(vertex_new(:,1)./(1+vertex_new(:,3)), vertex_new(:,2)./(1+vertex_new(:,3)));
[~, index] = sort(abs(z(face(:,1)))+abs(z(face(:,2)))+abs(z(face(:,3))));
ni = index(2); % since index(1) must be bi, not really
ns = sum(abs(z(face(bi,[1 2 3]))-z(face(bi,[2 3 1]))))/3;
ss = sum(abs(w(face(ni,[1 2 3]))-w(face(ni,[2 3 1]))))/3;
z = z*(sqrt(ns*ss))/ns;
dz2 = abs(z).^2;
vertex_new = [2*real(z)./(1+dz2), 2*imag(z)./(1+dz2), (-1+dz2)./(1+dz2)];

% south pole stereographic projection
uv = [vertex_new(:,1)./(1+vertex_new(:,3)),vertex_new(:,2)./(1+vertex_new(:,3))];
mu = compute_bc(face,uv,vertex);
% find the south pole
[~,ind] = sort(vertex_new(:,3));
nv = size(vertex,1);
fixed = ind(1:min(nv,50));
% reconstruct map with given mu and some fixed point
[fuv,fmu] = linear_beltrami_solver(face,uv,mu,fixed,uv(fixed,:));
fz = complex(fuv(:,1),fuv(:,2));
dfz2 = abs(fz).^2;
uvw = [2*real(fz)./(1+dfz2), 2*imag(fz)./(1+dfz2), -(-1+dfz2)./(1+dfz2)];