Matlab: Algorithm for voronoi diagram of ellipses

Question

Are there any algorithms to implement Voronoi diagram that bounds the ellipses? The diagram would look like the pictures here voronoi diagram of ellipses

Can anyone share some links,tutorials,codes etc related to it?

Thanks in advance.

Rachal

Answer 1

Here's an algorithm that uses the distance transform together with the watershed algorithm to draw a Voronoi diagram for ellipses.

%# first, define some ellipses (for simplicity, I use 0/90 orientation)
ellipses = [10,20,5,10;30,10,10,7;40,40,8,3];

%# put the ellipses into an image (few pixels, therefore pixelated)
img = false(50);
[xx,yy]=ndgrid(1:50,1:50);
for e = 1:size(ellipses,1),img = img | (xx-ellipses(e,1)).^2/ellipses(e,3)^2 + (yy-ellipses(e,2)).^2/ellipses(e,4)^2 <= 1;end

%# perform the distance transform
dt = bwdist(img);

%# apply the watershed algorithm. 
%# ws==0 are the lines for the Voronoi diagram
ws = watershed(dt);

%# create a RGB image and display
%# note: for yellow lines, replace the last
%# "ws==0" by "zeros(size(ws))", so that you
%# only put ws into the red and green channel (=yellow)
rgb = cat(3,ws==0,ws==0,ws==0)); 
%# add the ellipses into the red channel
rgb(:,:,1) = rgb(:,:,1) | img;
imshow(rgb)

Merissameristem

What can i do to the codes to control the color of line and ellipses? – Rachal

@Ivy: rgb = repmat(ws==0,1,1,3);rgb(:,:,1) = rgb(:,:,1) | img;imshow(rgb) – Merissameristem

How can i show both ellipses and lines at once in image? Because i imshow(dt), only ellipses display, if i imshow(ws), only lines. and i tried the rgb, there was an error for 'repmat': too many input arguments. What could i do? – Rachal

@Ivy: sorry, typo. It should be rgb=repmat(ws==0,[1 1 3]);. – Merissameristem

:sorry, i still have no idea how to display the image, using'imshow(???)'. Assume if background is black, ellipse' color is red, line color is yellow. Could you please add some codes in detail? Thanks a lot. – Rachal

Answer 2

Just in case, this is an example from the Mathematica help system:

(*Generate ellipses*)
p= Rasterize@Graphics@Table[
          Rotate[
              Disk[RandomReal[10, 2],          (*Rnd position*)
                   RandomReal[{.3, 1.5}, 2]],  (*Rnd radii*)
          RandomReal[Pi]], {i, 10}]            (*Rnd rotation*)

(*Compute Voronoi*)

LaplacianGaussianFilter[DistanceTransform[p], 2] // ImageAdjust

It is not an exact calculation, but fair enough for practical uses.

Crooked

Answer 3

Based on your recent track of questions, I understand that you've been working on drawing rasterized ellipses on top of an RGB image. You would like to be able to specify the ellipses location, shape and color. You wanted the ellipses to be clipped at the boundaries, and also to be non-overlapping. Now you are looking to draw the lines that divides the space in a manner similar to Voronoi diagrams (but with ellipses instead of points).

For this particular question, as @Jonas showed, the solution is to use the distance transform together with the watershed algorithm.

I thought I continue with my previous example, and extend it with Jonas's idea, to showcase the entire process. Hope you find it useful..

The code uses the calculateEllipse function to compute the coordinates of the points that make up an ellipse, as well as imoverlay function for setting specified pixels of an image to some chosen color.

%# color image (canvas to draw on)
I = imread('pears.png');
sz = size(I);

%# random ellipses
num = 20;
centers = bsxfun(@times, rand(num,2), sz([2 1]));   %# center x/y-coords
radii = bsxfun(@times, rand(num,2), [300 50])+10;   %# major/minor axis length
angles = rand(num,1) .* 360;                        %# angle of rotation
ex = cell(num,1);                                   %# vertices x-coords
ey = cell(num,1);                                   %# vertices y-coords

%# label image, used to hold rasterized ellipses
L = zeros(sz(1),sz(2));

%# randomly place ellipses one-at-a-time, skip if overlaps previous ones
flag = false(num,1);
for i=1:num
    %# ellipse we would like to draw directly on image matrix
    [ex{i},ey{i}] = calculateEllipse(centers(i,1),centers(i,2), ...
        radii(i,1),radii(i,2), angles(i), 100);

    %# create mask for image pixels inside the ellipse polygon
    mask = poly2mask(ex{i},ey{i}, sz(1),sz(2));

    %# check if there is no existing overlapping ellipse
    if all( L(mask)==0 )
        %# use the mask to place the ellipse in the label image
        L(mask) = sum(flag)+1;    %# assign value using an increasing counter
        flag(i) = true;
    end
end

%# filter ellipses to only those that made through the overlap test
num = sum(flag);
centers = centers(flag,:);
radii = radii(flag,:);
angles = angles(flag);
ex = ex(flag);
ey = ey(flag);

%# rasterized voroni diagram of the ellipses [Jonas]
E = (L ~= 0);                             %# ellipses as binary image
WS = watershed( bwdist(E) );              %# distance transform + watershed
WS = (WS == 0);                           %# WS==0 corresponds voronoi diagram
WS = bwmorph(WS, 'thicken',1);            %# thicken the lines

%# set pixels corresponding to voronoi diagram to white
II = I;
II = imoverlay(II, WS, [1 1 1]);          %# you can customize the color here

%# set pixels corresponding to ellipses using specified colors
clr = hsv(num);                           %# color of each ellipse
for i=1:num
    mask = bwperim(L==i,8);               %# get perimeter of the ellipse mask
    mask = bwmorph(mask, 'thicken',1);    %# thicken the ellipse perimeter
    II = imoverlay(II, mask, clr(i,:));   %# set those pixels with RGB color
end

%# show final rasterized image (image + ellipses + voronoi diagram)
figure, imshow(II, 'InitialMagnification',100, 'Border','tight')

Jephum

I tried to draw the solid ellipse, so i took out the line 'mask = bwperim(L==i,8);'which is to get the perimeter of ellipse, but the output came out with one solid ellipse only.Why did this happen?Is it some kind of logical error? – Rachal

@Ivy: if you want to draw solid ellipses, you should also remove the line immediately after it (thicken the perimeter), thus you would only have: mask = (L==i); – Jephum

Yes, after remove the two lines, change 'II = imoverlay(II, mask, clr(i,:));' to ' II = imoverlay(II, L, clr(i,:));' will work also. By the way, if my I is 'I=zeros(500,500,3)', how to change to background color from black to (assume) cyan[1 0 1]? Should i need to apply the imoverlay again? or there is a fast way? – Rachal

@Ivy: To be exact, it should be: II = imoverlay(II, L==i, clr(i,:));. Now if you were starting with a black image I(:) = 0;, you can use I = imoverlay(I,true(sz(1),sz(2)),[1 0 1]); to change its background color to cyan. – Jephum

I can't wrap my head around the part for checking overlapping ellipses,especially this line if all( L(mask)==0 ). what does it mean? – Rachal

@Ivy: basically in the step just before, we build a mask containing a logical matrix (true/false) that indicates the inside of the ellipse. Now before we place it in the label matrix, we check it doesn't overlap any existing ellipses, this is done by checking if all of the "true" pixels indicated by the mask all have zero values in the label matrix (if at least one pixel was non-zero, it means we already placed an ellipse passing by that location, and thus we skip the current one) – Jephum

:Thks.What if the ellipses just touch, how can i check it and skip the latter one also? Because when i increase the num to around 50, while i set all radius at the same value [10 6], touching happens during several trivals.And after bwdist and watershed,the two touching ellipses will be treated as one and thus get enclosed(two ellipses in a cell).That's not what i expect. – Rachal

@Ivy: ah I see your problem.. One possibility is to thicken the ellipse by 1 just before doing the overlap-test (a copy of it to be exact), that way we are additionally testing a perimeter one-pixel thick around its neighborhood, thus avoiding the "barely touching" case. So make the following changes: mask2 = bwmorph(mask, 'thicken',1); then do the test with this changed mask instead: if all( L(mask2)==0 ) ... – Jephum

@Ivy: compare the result on this example: before and after – Jephum

Answer 4

I dont know what you mean with "ellipses" here. But there is a implementation for voronoi diagrams of in C++ by Stephan Fortune / Shane O'Sullivan,

http://www.skynet.ie/~sos/mapviewer/voronoi.php

Iong

Ellipses take place of the 'points'(x,y) that are bounded by the lines. – Rachal