Code covered by the BSD License  

Highlights from
3D Voronoi diagram

image thumbnail

3D Voronoi diagram

by

 

13 May 2013 (Updated )

Given a list of centroids (.txt), it computes the 2D/3D voronoi diagram.

VoronoiCells.m
% Creates 3-D voronoi image given a centroid list (.txt)
% and extracts each Voronoi cell from the original image

% list of centroid locations
%i.e.
%Centroids.txt
%   x          y      z
%---------------------------
%| 881.18	920.484	14.89   |
%| 729.081	882.986	29.403  |
%| 611.288	322.088	25.509  |
%| 574.717	859.578	14.372  |
%| 64.306	743.749	33.244  |
%| 316.106	346.42	30.635  |
%---------------------------

% For 2D implementation, set z column to 1

%Set up
centroids = load('Centroids.txt');
originalImage='montage_kt01710_w311GFPdsu_BS_CV-_2000_3500_1000_1000_50-170_crop.tif';
    
%% Label Image
labelFilename = 'labelImage.tif';
InfoImage=imfinfo(originalImage);
mImage=InfoImage(1).Width;
nImage=InfoImage(1).Height;
NumberImages=length(InfoImage);
intensityImage=zeros(nImage,mImage,NumberImages,'uint8');
 
for i=1:NumberImages
   intensityImage(:,:,i)=imread(originalImage,'Index',i,'Info',InfoImage);
end

% check if label (voronoi) image already exist else create label image
if exist(labelFilename,'file') == 2
    disp('Loading existing label image...')
    InfoImage=imfinfo(labelFilename);
    mImage=InfoImage(1).Width;
    nImage=InfoImage(1).Height;
    NumberImages=length(InfoImage);
    labelImage=zeros(nImage,mImage,NumberImages,'uint8');

    for i=1:NumberImages
       labelImage(:,:,i)=imread(labelFilename,'Index',i,'Info',InfoImage);
    end
else

    labelImage = zeros(size(intensityImage));

    %find closest centroid for each voxel
    eucDistance = zeros(size(centroids,1),1);
    for i = 1:size(intensityImage,1)
        for j = 1:size(intensityImage,2)
            for k = 1:size(intensityImage,3)
                cur_idx = [i j k];
                eucDistance = sqrt((centroids(:,1)-i).^2+(centroids(:,2)-j).^2+(centroids(:,3)-k).^2);
                [~,index] = min(eucDistance);
                labelImage(i,j,k)=uint8(index);
            end
        end
    end

    for k = 1:size(intensityImage,3)
        imwrite(uint8(labelImage(:,:,k)),labelFilename,'WriteMode','append','Compression','none');
    end
    disp('Voronoi image Saved')

end

%% Extract each Voronoi cell
for cell_num = 1:size(centroids,1)
    [x y z] = ind2sub(size(labelImage),find(labelImage == cell_num));
    min_x = min(x);    max_x = max(x);
    min_y = min(y);    max_y = max(y);
    min_z = min(z);    max_z = max(z);

    cropImage = intensityImage(min_x:max_x,min_y:max_y,min_z:max_z);
    
    cellFilename = ['cropImage_',num2str(cell_num),'.tif'];
    for k = 1:size(cropImage,3)
        imwrite(uint8(cropImage(:,:,k)),cellFilename,'tif','WriteMode','append','Compression','none');
    end
    disp([cellFilename,' saved'])
end
disp('Extraction done')

Contact us