Tri mesh within a voronoi

I have a rectangular boundary with an internal voronoi diagram. Does anyone have any ideas on an alternate method for developing a tri mesh within the cells. I have been working on one with little avail. I start with the vx, and vy output from the voronoi (plus the boundary nodes for all voronoi boundary intersections). e.g.

[vx,vy] = voronoi(dt);

From here I attempt to setup to use the mesh2d toolbox from file exchange (identifying the nodes, edges and faces). My problem is my method keeps missing some of the cells and I cannot track down the problem. I am open to new code suggestions, but I have included my code just in case.

function [panel_nodes panel_mesh_map node edge] = Create_mesh(all_vx,all_vy,x1,x2,y1,y2,options,hdata)

% Mesh Cells

% Collect all nodes and remove duplicates
node=zeros(2*length(all_vx),2);
edge=zeros(length(all_vx),2);
for j = 1:length(all_vx)
    node(2*j-1,1) = all_vx(1,j);
    node(2*j-1,2) = all_vy(1,j);
    node(2*j,1) = all_vx(2,j);
    node(2*j,2) = all_vy(2,j);
    edge(j,1) = 2*j-1;
    edge(j,2) = 2*j;
end

num_node = length(node);
% Add the boundary nodes
node(num_node+1,1) = x1;
node(num_node+1,2) = y1;
node(num_node+2,1) = x2;
node(num_node+2,2) = y1;
node(num_node+3,1) = x2;
node(num_node+3,2) = y2;
node(num_node+4,1) = x1;
node(num_node+4,2) = y2;

% Remove duplicate nodes
[node junk node_map] = unique(node,'rows');  

% node_map is the maping of original location(cell number) to the new location
% of the node(value of the cell)
for j = 1:length(edge)
    if edge(j,1)==edge(j,2)
        edge(j,:)=[];
    end
    edge(j,1) = node_map(edge(j,1));
    edge(j,2) = node_map(edge(j,2));
end

node = cell2mat(num2cell(node));

% Remove zero-length edges
changing_length = length(edge); % the length of edge will change as rows are removed
k = 0;
for j = 1:length(edge)
    k = k+1;
    if k<changing_length
    if edge(k,1)==edge(k,2)
        edge(k,:)=[];
        k = k-1; % next row now equal to j so check must occur for same row again
        changing_length = changing_length-1;
    end
    end
end

% Add the boundary edges

% setup indicies of the Y components in order for edges
y1_indx = 1;
y2_indx = 1;
for j = 1:length(node)
    if node(j,2) == y1
        y1_stor(y1_indx) = j;
        y1_indx = y1_indx + 1;
    elseif node(j,2) == y2
        y2_stor(y2_indx) = j;
        y2_indx = y2_indx + 1;
    end
end

for j = 1:length(y1_stor)-1
    edge(length(edge)+1,1) = y1_stor(j);
    edge(length(edge),2) = y1_stor(j+1);
end

for j = 1:length(y2_stor)-1
    edge(length(edge)+1,1) = y2_stor(j);
    edge(length(edge),2) = y2_stor(j+1);
end

for j = 1:length(node)
    if node(j,1) == x1 && node(j,2) < y2 % left side
        edge(length(edge)+1,1) = j;
        edge(length(edge),2) = j+1;
    elseif node(j,1) == x2 && node(j,2) < y2 % right side
        edge(length(edge)+1,1) = j;
        edge(length(edge),2) = j+1;
    end
end

% Setup the faces of the cells
% Calculate all counter-clockwise
face = [];
for j = 1:length(edge);
    shared_edge = [0 0 0 0];
    face_stor = j;
    edge_stor = [];
    indx = 1;
    initial = 1;
    while shared_edge(indx,1) ~= edge(j,1)
        if initial == 1
        new_edge(1,1) = edge(j,1);
        new_edge(1,2) = edge(j,2);
        elseif shared_edge(indx,4) == 0 %no change needed
        new_edge(1,1) = shared_edge(indx,1);
        new_edge(1,2) = shared_edge(indx,2);
        elseif shared_edge(indx,4) == 1 %needs to change order
        new_edge(1,1) = shared_edge(indx,2);
        new_edge(1,2) = shared_edge(indx,1);
        end
        shared_edge = [];
        initial = 0;
    for k = 1:length(edge); % find edges that share node 2 of primary edge and orient it to point to node 2
      if k~=j
        if edge(k,1) == new_edge(1,2)
            shared_edge = [shared_edge;edge(k,2),edge(k,1),k,1]; % [node1, node 2, row in edge, (1or0) orientation]
        elseif edge(k,2) == new_edge(1,2)
            shared_edge = [shared_edge;edge(k,1),edge(k,2),k,1];
        end
      end
    end

%     Once all shared edges are found and oriented, calculate angles to
%     find smallest positive angle, counter-clock rotation
    for k = 1:length(shared_edge(:,1))
%       ang = mod( atan2( (x2-x1)*(y4-y3)-(y2-y1)*(x4-x3) ,(x2-x1)*(x4-x3)+(y2-y1)*(y4-y3) ) , 2*pi)
      angle = mod( atan2( (node(new_edge(1,2),1)-node(new_edge(1,1),1))*(node(shared_edge(k,2),2)-node(shared_edge(k,1),2))-(node(new_edge(1,2),2)-node(new_edge(1,1),2))*(node(shared_edge(k,2),1)-node(shared_edge(k,1),1)),(node(new_edge(1,2),1)-node(new_edge(1,1),1))*(node(shared_edge(k,2),1)-node(shared_edge(k,1),1))+(node(new_edge(1,2),2)-node(new_edge(1,1),2))*(node(shared_edge(k,2),2)-node(shared_edge(k,1),2)) ) , 2*pi); 
%        if angle == 0
%            angle = 0.0001;
%        end
      shared_edge(k,5) = angle;
    end
    shared_edge_ang = shared_edge(:,5)';
    max_ang = max(shared_edge_ang((shared_edge_ang>0)));
    indx = min(find(shared_edge_ang==max_ang));
% Store edge as part of closed face - End if connected with primary edge
    if shared_edge(indx,5) >= 3.141592653589793 && shared_edge(indx,5) <= 3.1415926535897932 % circling outer boundary
        face_stor = [];
        break
    end
    face_stor = [face_stor,shared_edge(indx,3)];
    edge_stor = [edge_stor,shared_edge(indx,2),shared_edge(indx,1)];
    end
%     Collect all complete faces
%     face{j+length(edge)} = face_stor;
    face{j} = face_stor;
end

% Remove similar faces
n=0;
changing_length = length(face);
for j = 1:length(face)
  n = n + 1;
  if n <= changing_length
    face_sorted_n = sort(face{n});
    m = j;
    for k = j+1:length(face)
        m = m + 1;
        changing_length = length(face);
        if m <= changing_length
          face_sorted_m = sort(face{m});
          if length(face_sorted_m) == length(face_sorted_n)
           if face_sorted_m == face_sorted_n
            face(m) = [];
            changing_length = changing_length - 1;
            m = m - 1;  % cells shifted so need to check same value again
           end
          end
        end
    end
  end
end

% Remove empty faces
k=0;
changing_length = length(face);
for j = 1:length(face)
    k = k +1;
    if k <= changing_length
      if isempty(face{k})
        face(k) = [];
        k = k-1; % cells shifted so need to check same value again
        changing_length = changing_length - 1;
      end
    end
end

% Generate Mesh

[panel_nodes panel_mesh_map] = meshfaces(node,edge,face,hdata,options);