function [ fout,vout ] = subdividemesh( f,v )
%[fout,vout] SUBDIVIDEMESH(f,v) subdivides the mesh once.
%   Detailed explanation goes here

%% Find all the unique edges

e1 = f(:,[1,2]);
e2 = f(:,[1,3]);
e3 = f(:,[2,3]);
edges = [e1
    e2
    e3];
edges = [min(edges(:,1),edges(:,2)) max(edges(:,1),edges(:,2))];
[edges,foo,edge_numbers] = unique(edges,'rows');


%% Find all the mdpts
midpoints = 0.5*(v(edges(:,1),:)+ v(edges(:,2),:));

vout = [v
    midpoints];


%% Get mdpt labels
n=size(f,1);
m=size(v,1);

m1=edge_numbers(1:n,1)+m;
m2=edge_numbers(n+1:2*n,1)+m;
m3=edge_numbers(2*n+1:3*n,1)+m;


f1=[f(:,1) m1 m2];
f2=[f(:,2) m1 m3];
f3=[f(:,3) m2 m3];
f4=[m1 m2 m3];
fout=[f1
      f2
      f3
      f4];