%function returs stlPoints fromat and ABC format if its needed,if not - just delete it and adopt to your needs
function [stlPoints Apoints Bpoints Cpoints ]= sphereTriangulation(numIterations, radius)
if nargin <1
radius = 1; % unit circle
numIterations = 3; % enough for testing
elseif nargin<2 %set default
radius = 1;
end%no more error handlers from me since I am not supposed to know for what You will use it :)
%basic Octahedron reff:http://en.wikipedia.org/wiki/Octahedron
% ( ?1, 0, 0 )
% ( 0, ?1, 0 )
% ( 0, 0, ?1 )
A=[1 0 0]*radius;
B=[0 1 0]*radius;
C=[0 0 1]*radius;
%from +-ABC crete initial triangles which define oxahedron
triangles=[ A; B; C;...
A; B; -C;....
% -x, +y, +-Z quadrant
-A; B; C;...
-A; B; -C;...
% -x, -y, +-Z quadrant
-A; -B; C;...
-A; -B; -C;...
% +x, -y, +-Z quadrant
A; -B; C;...
A; -B; -C];%-----STL-similar format
%for simplicity lets break into ABC points...
selector = 1:3:numel(triangles(:,1))-1 ;
Apoints = triangles(selector ,:) ;
Bpoints = triangles(selector+1,:) ;
Cpoints = triangles(selector+2,:) ;
%in every of numIterations
for iteration = 1:numIterations
%devide every of triangle on three new
% ^ C
% / \
% AC/2 /_4_\CB/2
% /\ 3 /\
% / 1\ /2 \
% A /____V____\B 1st 2nd 3rd 4th
% AB/2
%new triangleSteck is [ A AB/2 AC/2; AB/2 B CB/2; AC/2 AB/2 CB/2 AC/2 CB/2 C]
AB_2 = (Apoints+Bpoints)/2;
%do normalization of vector
AB_2 = arsUnit(AB_2,radius);...same for next 2 lines
AC_2 = (Apoints+Cpoints)/2; AC_2 = arsUnit(AC_2,radius);
CB_2 = (Cpoints+Bpoints)/2; CB_2 = arsUnit(CB_2,radius);
Apoints = [ Apoints;... A point from 1st triangle
AB_2;... A point from 2nd triangle
AC_2;... A point from 3rd triangle
AC_2];... A point from 4th triangle..same for B and C
Bpoints = [AB_2; Bpoints; AB_2; CB_2 ];
Cpoints = [AC_2; CB_2 ; CB_2; Cpoints];
end
%now tur points back to STL-like format....
numPoints = numel(Apoints(:,1)) ;
selector = 1:numPoints ;
selector = selector(:) ;
selector = [selector, selector+numPoints, selector+2*numPoints];
selector = selector' ;
selector = selector(:) ;%thats it
stlPoints = [Apoints; Bpoints; Cpoints] ;
stlPoints = stlPoints(selector,:) ;
end
%vectorized norm() function
function [rez]=arsNorm(A)
rez = A(:,1).*A(:,1) + A(:,2).*A(:,2) + A(:,3).*A(:,3);
rez = sqrt(rez);
end
%vectorized unit() functon
function [rez]=arsUnit(A,radius)
normOfA = arsNorm(A) ;
rez = A ./ [normOfA normOfA normOfA];
rez = rez * radius ;
end