sphere_tri - generate a triangle mesh approximating a sphere
Usage: FV = sphere_tri(shape,Nrecurse,r,winding)
shape is a string, either of the following:
'ico' starts with icosahedron (most even, default)
'oct' starts with octahedron
'tetra' starts with tetrahedron (least even)
Nrecurse is int >= 0, setting the recursions (default 0)
r is the radius of the sphere (default 1)
winding is 0 for clockwise, 1 for counterclockwise (default 0). The
matlab patch command gives outward surface normals for clockwise
order of vertices in the faces (viewed from outside the surface).
FV has fields FV.vertices and FV.faces. The vertices
are listed in clockwise order in FV.faces, as viewed
from the outside in a RHS coordinate system.
The function uses recursive subdivision. The first
approximation is an platonic solid, either an icosahedron,
octahedron or a tetrahedron. Each level of refinement
subdivides each triangle face by a factor of 4 (see also
mesh_refine). At each refinement, the vertices are
projected to the sphere surface (see sphere_project).
A recursion level of 3 or 4 is a good sphere surface, if
gouraud shading is used for rendering.
The returned struct can be used in the patch command, eg:
% create and plot, vertices: [2562x3] and faces: [5120x3]
FV = sphere_tri('ico',4,1);
lighting phong; shading interp; figure;
patch('vertices',FV.vertices,'faces',FV.faces,...
'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]);
axis off; camlight infinite; camproj('perspective');
See also: mesh_refine, sphere_project0001 function [FV] = sphere_tri(shape,maxlevel,r,winding) 0002 0003 % sphere_tri - generate a triangle mesh approximating a sphere 0004 % 0005 % Usage: FV = sphere_tri(shape,Nrecurse,r,winding) 0006 % 0007 % shape is a string, either of the following: 0008 % 'ico' starts with icosahedron (most even, default) 0009 % 'oct' starts with octahedron 0010 % 'tetra' starts with tetrahedron (least even) 0011 % 0012 % Nrecurse is int >= 0, setting the recursions (default 0) 0013 % 0014 % r is the radius of the sphere (default 1) 0015 % 0016 % winding is 0 for clockwise, 1 for counterclockwise (default 0). The 0017 % matlab patch command gives outward surface normals for clockwise 0018 % order of vertices in the faces (viewed from outside the surface). 0019 % 0020 % FV has fields FV.vertices and FV.faces. The vertices 0021 % are listed in clockwise order in FV.faces, as viewed 0022 % from the outside in a RHS coordinate system. 0023 % 0024 % The function uses recursive subdivision. The first 0025 % approximation is an platonic solid, either an icosahedron, 0026 % octahedron or a tetrahedron. Each level of refinement 0027 % subdivides each triangle face by a factor of 4 (see also 0028 % mesh_refine). At each refinement, the vertices are 0029 % projected to the sphere surface (see sphere_project). 0030 % 0031 % A recursion level of 3 or 4 is a good sphere surface, if 0032 % gouraud shading is used for rendering. 0033 % 0034 % The returned struct can be used in the patch command, eg: 0035 % 0036 % % create and plot, vertices: [2562x3] and faces: [5120x3] 0037 % FV = sphere_tri('ico',4,1); 0038 % lighting phong; shading interp; figure; 0039 % patch('vertices',FV.vertices,'faces',FV.faces,... 0040 % 'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]); 0041 % axis off; camlight infinite; camproj('perspective'); 0042 % 0043 % See also: mesh_refine, sphere_project 0044 % 0045 0046 0047 0048 % $Revision: 1.2 $ $Date: 2005/07/20 23:07:03 $ 0049 0050 % Licence: GNU GPL, no implied or express warranties 0051 % Jon Leech (leech @ cs.unc.edu) 3/24/89 0052 % icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93 0053 % 06/2002, adapted from c to matlab by Darren.Weber_at_radiology.ucsf.edu 0054 % 05/2004, reorder of the faces for the 'ico' surface so they are indeed 0055 % clockwise! Now the surface normals are directed outward. Also reset the 0056 % default recursions to zero, so we can get out just the platonic solids. 0057 % 0058 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0059 0060 eegversion = '$Revision: 1.2 $'; 0061 fprintf('SPHERE_TRI [v %s]\n',eegversion(11:15)); tic 0062 0063 if ~exist('shape','var') || isempty(shape), 0064 shape = 'ico'; 0065 end 0066 fprintf('...creating sphere tesselation based on %s\n',shape); 0067 0068 % default maximum subdivision level 0069 if ~exist('maxlevel','var') || isempty(maxlevel) || maxlevel < 0, 0070 maxlevel = 0; 0071 end 0072 0073 % default radius 0074 if ~exist('r','var') || isempty(r), 0075 r = 1; 0076 end 0077 0078 if ~exist('winding','var') || isempty(winding), 0079 winding = 0; 0080 end 0081 0082 0083 % ----------------- 0084 % define the starting shapes 0085 0086 shape = lower(shape); 0087 0088 switch shape, 0089 case 'tetra', 0090 0091 % Vertices of a tetrahedron 0092 sqrt_3 = 0.5773502692; 0093 0094 tetra.v = [ sqrt_3, sqrt_3, sqrt_3 ; % +X, +Y, +Z - PPP 0095 -sqrt_3, -sqrt_3, sqrt_3 ; % -X, -Y, +Z - MMP 0096 -sqrt_3, sqrt_3, -sqrt_3 ; % -X, +Y, -Z - MPM 0097 sqrt_3, -sqrt_3, -sqrt_3 ]; % +X, -Y, -Z - PMM 0098 0099 % Structure describing a tetrahedron 0100 tetra.f = [ 1, 2, 3; 0101 1, 4, 2; 0102 3, 2, 4; 0103 4, 1, 3 ]; 0104 0105 FV.vertices = tetra.v; 0106 FV.faces = tetra.f; 0107 0108 case 'oct', 0109 0110 % Six equidistant points lying on the unit sphere 0111 oct.v = [ 1, 0, 0 ; % X 0112 -1, 0, 0 ; % -X 0113 0, 1, 0 ; % Y 0114 0, -1, 0 ; % -Y 0115 0, 0, 1 ; % Z 0116 0, 0, -1 ]; % -Z 0117 0118 % Join vertices to create a unit octahedron 0119 oct.f = [ 1 5 3 ; % X Z Y - First the top half 0120 3 5 2 ; % Y Z -X 0121 2 5 4 ; % -X Z -Y 0122 4 5 1 ; % -Y Z X 0123 1 3 6 ; % X Y -Z - Now the bottom half 0124 3 2 6 ; % Y Z -Z 0125 2 4 6 ; % -X Z -Z 0126 4 1 6 ]; % -Y Z -Z 0127 0128 FV.vertices = oct.v; 0129 FV.faces = oct.f; 0130 0131 case 'ico', 0132 0133 % Twelve vertices of icosahedron on unit sphere 0134 tau = 0.8506508084; % t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2) 0135 one = 0.5257311121; % one=1/sqrt(1+t^2) , unit sphere 0136 0137 ico.v( 1,:) = [ tau, one, 0 ]; % ZA 0138 ico.v( 2,:) = [ -tau, one, 0 ]; % ZB 0139 ico.v( 3,:) = [ -tau, -one, 0 ]; % ZC 0140 ico.v( 4,:) = [ tau, -one, 0 ]; % ZD 0141 ico.v( 5,:) = [ one, 0 , tau ]; % YA 0142 ico.v( 6,:) = [ one, 0 , -tau ]; % YB 0143 ico.v( 7,:) = [ -one, 0 , -tau ]; % YC 0144 ico.v( 8,:) = [ -one, 0 , tau ]; % YD 0145 ico.v( 9,:) = [ 0 , tau, one ]; % XA 0146 ico.v(10,:) = [ 0 , -tau, one ]; % XB 0147 ico.v(11,:) = [ 0 , -tau, -one ]; % XC 0148 ico.v(12,:) = [ 0 , tau, -one ]; % XD 0149 0150 % Structure for unit icosahedron 0151 ico.f = [ 5, 8, 9 ; 0152 5, 10, 8 ; 0153 6, 12, 7 ; 0154 6, 7, 11 ; 0155 1, 4, 5 ; 0156 1, 6, 4 ; 0157 3, 2, 8 ; 0158 3, 7, 2 ; 0159 9, 12, 1 ; 0160 9, 2, 12 ; 0161 10, 4, 11 ; 0162 10, 11, 3 ; 0163 9, 1, 5 ; 0164 12, 6, 1 ; 0165 5, 4, 10 ; 0166 6, 11, 4 ; 0167 8, 2, 9 ; 0168 7, 12, 2 ; 0169 8, 10, 3 ; 0170 7, 3, 11 ]; 0171 0172 FV.vertices = ico.v; 0173 FV.faces = ico.f; 0174 end 0175 0176 0177 % ----------------- 0178 % refine the starting shapes with subdivisions 0179 if maxlevel, 0180 0181 % Subdivide each starting triangle (maxlevel) times 0182 for level = 1:maxlevel, 0183 0184 % Subdivide each triangle and normalize the new points thus 0185 % generated to lie on the surface of a sphere radius r. 0186 FV = mesh_refine_tri4(FV); 0187 FV.vertices = sphere_project(FV.vertices,r); 0188 0189 % An alternative might be to define a min distance 0190 % between vertices and recurse or use fminsearch 0191 0192 end 0193 end 0194 0195 if winding, 0196 fprintf('...returning counterclockwise vertex order (viewed from outside)\n'); 0197 FV.faces = FV.faces(:,[1 3 2]); 0198 else 0199 fprintf('...returning clockwise vertex order (viewed from outside)\n'); 0200 end 0201 0202 t=toc; fprintf('...done (%6.2f sec)\n\n',t); 0203 0204 return