How to generate Voronoi diagram using Quasi-Euclidean distance transform?

Dear colleagues, I am trying to create Voronoi diagrams using Quasi-Euclidean distance transform (QEDT), which is easily obtained in matlab using the bwdist(.) command. I got some codes from http://math.stackexchange.com/questions/960366/how-do-i-compose-a-voronoi-diagram that helped me to get started. The problem is that I am not confident with the result I get. I thought that QEDT could give a more accurate solution, just close enough to the ideal Euclidean distance. Below is the code I have edited, and I don't know if it correct--provided a wide deviation between the two distances. Please check the comment lines for details. My edits are from line 13 to 29.

close all;
% Voronoi script
 p=[100,20;100,190;10,100;170,100;100,100;155,110];
voronoi(p(:,1),p(:,2))
axis([1,200,1,200])
n=zeros (200);
m=n;
d=zeros(size(p));
for i=1:200
    for j=1:200
        % euclidean
        d1=(sum((p-repmat([i,j],size(p,1),1)).^2,2)).^(0.5);   
        % Quasi-Euclidean
        % The formula for the Quasi-Euclidean distance is given as:
        % abs(x1-x2) + (sqrt(2)-1)*abs(y1-y2),  if abs(x1-x2) > abs(y1-y2)
        % (sqrt(2)-1)*abs(x1-x2) + abs(y1-y2),  otherwise
        % In matlab, typing bwdist in the command window gives details
        % about several distances transformations, including
        % Quasi-Euclidean
        q = repmat([i,j],size(p,1),1);
        x1 = p(:,1);
        x2 = q(:,1);
        y1 = p(:,2);
        y2 = q(:,2);
        if (abs(x1 - x2) > abs(y1 - y2))
            d2 = sum((abs(x1 - x2) + (sqrt(2) - 1)*abs(y1 - y2)),2);
        else
            d2 = sum(((sqrt(2) - 1)*abs(x1 - x2) + abs(y1 - y2)),2);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Manhattan distance
        % d2 = sum(abs(p-repmat([i,j],size(p,1),1)),2);
        idx1=find(d1==min(d1));
        idx2=find(d2==min(d2));
        if length(idx1)>1
            n(i,j)=0;
        else
            n(i,j)=idx1(1);
        end
        if length(idx2)>1
            m(i,j)=0;
        else
            m(i,j)=idx2(1);
        end
    end
end
figure;
axis([1,200,1,200])
%imagesc(flip(n'));
imagesc(n');
hold on;
pf=p;
scatter(pf(:,1),pf(:,2));

figure;
axis([1,200,1,200])
imagesc(m');
hold on;
pf=p;
scatter(pf(:,1),pf(:,2));