function [points,edges,faces]=splitec(points,edges,faces)
    fnew=[];
    enew=[];
    ptnew=[];
    ptl=size(points,1);
    etl=size(edges,1);
    edges=[edges zeros(etl,2)];
    for fi=1:size(faces,1)
        fe=faces(fi,:);
        eu=zeros(6,2);
        pu=zeros(6,3);
        for ei=1:3
            if fe(ei)>0 //traverse edge forward
                neg=%f;
            else
                neg=%t;
                fe(ei)=-fe(ei);
            end
            ep=edges(fe(ei),:);
            if ep(3)<1 //edge has not been split yet
                ptm=[points(ep(1),:);points(ep(2),:)];
                ptm=mean(ptm,1);
                ptm=ptm/norm(ptm(1:3));
                //ptm(4)=ord;
                ptnew=[ptnew;ptm];
                a=ptl+size(ptnew,1);
                enew=[enew;ep(1) a;a ep(2)];
                eu(ei,:)=size(enew,1)+[-1 0];
                pu(ei,:)=[ep(1) a ep(2)];
                edges(fe(ei),3:4)=eu(ei,:);
            else
                eu(ei,:)=ep(3:4);
                pu(ei,:)=[enew(ep(3),1) enew(ep(3),2) enew(ep(4),2)];
            end
            if(neg)
                eu(ei,:)=[-eu(ei,2) -eu(ei,1)];
                pu(ei,:)=[pu(ei,3) pu(ei,2) pu(ei,1)];
            end
        end
        enew1=[pu(1,2) pu(3,2)];
        enew2=[pu(2,2) pu(1,2)];
        enew3=[pu(3,2) pu(2,2)];
        enew=[enew;enew1;enew2;enew3];
        eun=size(enew,1)+[-2 -1 0];
        fnew1=[eu(1,1) eun(1) eu(3,2)];
        fnew2=[eu(1,2) eu(2,1) eun(2)];
        fnew3=[eun(3) eu(2,2) eu(3,1)];
        fnew4=[eun(3) eun(2) eun(1)];
        fnew=[fnew;fnew1;fnew2;fnew3;fnew4];
    end
    points=[points; ptnew];
    edges=enew;
    faces=fnew;
endfunction
return
//spherical points distribution
//using the face-centered split method
//orthogonal axes octahedron
if %f
points=[
-1 0 0;
 1 0 0;
 0 -1 0;
 0 1 0;
 0 0 -1;
 0 0 1]; //orthogonal axes, order 0
edges=[
1 3 0 0;
1 4 0 0;
1 5 0 0;
1 6 0 0;
2 3 0 0;
2 4 0 0;
2 5 0 0;
2 6 0 0;
3 5 0 0;
3 6 0 0;
4 5 0 0;
4 6 0 0]; //lines as point indices, split 0
faces=[
 1 9 -3;
 1 10 -4;
 2 11 -3;
 2 12 -4;
 5 9 -7;
 5 10 -8;
 6 11 -7;
 6 12 -8]; //triangles as edge indices, negative = reverse edge
end
order_gen=1;
for ord=1:order_gen
    [points,edges,faces]=splitec(points,edges,faces);
end
//scatter3(points(:,1),points(:,2),points(:,3),[],points(:,4));
scatter3(points(:,1),points(:,2),points(:,3));
isoview
