function [points,edges,faces]=splitfc(points,edges,faces)
    fnew=[];
    enew=[];
    ptnew=[];
    ptl=size(points,1);
    etl=size(edges,1);
    for fi=1:size(faces,1)
        fe=faces(fi,:);
        pu=zeros(3,2);
        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),1:2);
            pu(ei,:)=ep;
            if(neg)
                pu(ei,:)=[pu(ei,2) pu(ei,1)];
            end
        end
        ptm=[points(pu(1,1),:);points(pu(2,1),:);points(pu(3,1),:)];
        mpt=mean(ptm,1);
        mpt=mpt/norm(mpt(1:3));
        //mpt(4)=ord;
        ptnew=[ptnew;mpt];
        a=ptl+fi;
        enew1=[pu(1,1) a];
        enew2=[pu(2,1) a];
        enew3=[pu(3,1) a];
        enew=[enew;enew1;enew2;enew3];
        eun=etl+size(enew,1)+[-2 -1 0];
        fnew1=[faces(fi,1) eun(2) -eun(1)];
        fnew2=[faces(fi,2) eun(3) -eun(2)];
        fnew3=[faces(fi,3) eun(1) -eun(3)];
        fnew=[fnew;fnew1;fnew2;fnew3];
    end
    points=[points; ptnew];
    edges=[edges; enew];
    faces=fnew;
endfunction
return
//spherical points distribution
//using the face-centered split method
if %t
points=[
-1 0 0 0;
 1 0 0 0;
 0 -1 0 0;
 0 1 0 0;
 0 0 -1 0;
 0 0 1 0]; //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]=splitfc(points,edges,faces);
end
scatter3(points(:,1),points(:,2),points(:,3));
isoview
