16 views (last 30 days)

Show older comments

Does anybody know a fast and accurate implementation for converting a binary image into a 2D triangulation? As an example consider the following image: http://tinypic.com/view.php?pic=25qulat&s=5. The code should be able to convert the left image into the right image...

I made an implementation myself, but to be honest, it is an (ugly) workaround which I prefer not to use anymore. However, it gets the job done in a small amount of time.

Here's an example on how my code works:

% generate binary image

nx = 100;

ny = 100;

image_binary = phantom('Modified Shepp-Logan', nx)>0;

% specify image domain

x = linspace(-1,1,nx);

y = -linspace(-1,1,ny);

% pad image with zeros in order to enable border at image boundaries

temp = zeros(size(image_binary)+2);

temp(2:end-1,2:end-1) = image_binary;

image_binary = temp;

x = [x(1)-(x(2)-x(1)), x, x(end)+(x(2)-x(1))];

y = [y(1)-(y(2)-y(1)), y, y(end)+(y(2)-y(1))];

[X,Y] = meshgrid(x,y);

% generate edge of the image (subtract eroded image from original image)

image_binary_edge = image_binary-imerode(image_binary,strel('disk',1));

% remove pixels with only one neighbour

image_binary_edge_filtered = imfilter(image_binary_edge,ones(3,3),'same');

image_binary_edge(image_binary_edge_filtered==2) = 0;

% calculate all connected components in image_binary_edge

cc = bwconncomp(image_binary_edge,8);

% initialize vectors for the delaunayTriangulation function

x_coor = [];

y_coor = [];

constraints = [];

max_dist = sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);

% loop over all components

for ii=1:cc.NumObjects

current = cc.PixelIdxList{ii};

x_coor_current = X(current);

y_coor_current = Y(current);

% reorder coordinates such that they are ordered in a clockwise fashion

x_coor_reordered = zeros(size(x_coor_current));

y_coor_reordered = zeros(size(y_coor_current));

x_coor_reordered(1) = x_coor_current(1);

y_coor_reordered(1) = y_coor_current(1);

x_coor_current(1) = [];

y_coor_current(1) = [];

kk=2;

while ~isempty(x_coor_current)

[index,dist] = knnsearch([x_coor_current,y_coor_current],[x_coor_reordered(kk-1),y_coor_reordered(kk-1)]);

% if dist is to large, than the current pixel is no neighbouring

% pixel, this is why we do not at these pixels to the reordered

% vectors

if(dist>2*max_dist)

x_coor_current(index) = [];

y_coor_current(index) = [];

else

x_coor_reordered(kk) = x_coor_current(index);

y_coor_reordered(kk) = y_coor_current(index);

x_coor_current(index) = [];

y_coor_current(index) = [];

kk = kk + 1;

end

end

x_coor_reordered = x_coor_reordered(1:kk-1); % remove zero entries

y_coor_reordered = y_coor_reordered(1:kk-1); % remove zero entries

% take only half of all border samples (this prevents oversampling of

% the border)

x_coor_reordered = x_coor_reordered(1:2:end);

y_coor_reordered = y_coor_reordered(1:2:end);

x_coor = [x_coor;x_coor_reordered];

y_coor = [y_coor;y_coor_reordered];

constraints_temp = [[length(constraints)+1:length(constraints)+length(x_coor_reordered)]',...

circshift([length(constraints)+1:length(constraints)+length(x_coor_reordered)]',-1)];

constraints = [constraints;constraints_temp];

end

% construct delaunay triangulation

dt = delaunayTriangulation(x_coor,y_coor,constraints);

% maintain only the interior

inside = dt.isInterior();

% Construct a triangulation that represents interior

tr = triangulation(dt(inside, :), dt.Points);

% at the moment, all vertices lie on the edge of the binary image,

% therefore, sample vertices inside the binary image as well:

pointstemp = tr.Points;

connectivityListtemp = tr.ConnectivityList;

pointsinside = zeros(size(X));

for t = 1:size(connectivityListtemp,1)

vertsXY = pointstemp(connectivityListtemp(t,:),:);

pointsinside = pointsinside | inpolygon(X,Y, vertsXY(:,1), vertsXY(:,2));

end

pointsinside(1:5:end,:) = 0;

pointsinside(2:5:end,:) = 0;

pointsinside(3:5:end,:) = 0;

pointsinside(4:5:end,:) = 0;

pointsinside(:,1:5:end) = 0;

pointsinside(:,2:5:end) = 0;

pointsinside(:,3:5:end) = 0;

pointsinside(:,4:5:end) = 0;

% construct the triangulation again

dt = delaunayTriangulation([x_coor;X(pointsinside==1)],[y_coor;Y(pointsinside==1)],constraints);

inside = dt.isInterior();

tr = triangulation(dt(inside, :), dt.Points);

% remove points which do not belong to triangle

Points = tr.Points;

ConnectivityList = tr.ConnectivityList;

ii=1;

while(ii<=length(Points))

if(~isempty(find(ConnectivityList == ii,1)))

ii = ii + 1;

else

Points(ii,:) = [];

ConnectivityList(ConnectivityList>ii) = ConnectivityList(ConnectivityList>ii)-1;

end

end

tr = triangulation(ConnectivityList,Points);

% plot the result

figure();

subplot(1,2,1)

imshow(image_binary,[])

title('Binary Image')

subplot(1,2,2)

triplot(tr.ConnectivityList,tr.Points(:,1),tr.Points(:,2))

title('triangulation')

Sven
on 28 Aug 2013

Edited: Sven
on 28 Aug 2013

Geert, here's how I'd do it. Note that I use isocontour for one step. Just a simple MATLAB "contour" call may also do the job, but that requires plotting to a figure so I went with an FEX function.

% Get a binary image

I = phantom('Modified Shepp-Logan', nx)>0;

% pad image with zeros in order to enable border at image boundaries

temp = zeros(size(I)+2);

temp(2:end-1,2:end-1) = I;

I = temp;

% Get an isocontour

contourThreshold = 0.5;

[Lines,Vertices,Objects] = isocontour(I,contourThreshold);

Vertices = fliplr(Vertices); % Get it back in XY from IJ

% Triangulate all pts in the isocontour and check which trias are in/out

DT = delaunayTriangulation(Vertices);

fc = DT.incenter;

in = interp2(I, fc(:,1), fc(:,2))>=contourThreshold;

% Show the result

figure,imagesc(I), hold on,

patch('vertices',DT.Points,'faces',DT.ConnectivityList(in,:),'FaceColor','g')

patch('vertices',DT.Points,'faces',DT.ConnectivityList(~in,:),'FaceColor','c')

plot(fc(in,1),fc(in,2),'b.', fc(~in,1),fc(~in,2),'y.')

for i=1:length(Objects)

Points=Objects{i};

plot(Vertices(Points,1),Vertices(Points,2),'Color','m');

end

Note that you could also get your vertices via bwperim rather than an isocontour... that would look like:

% Get an isocontour

[a,b] = find(bwperim(I));

Vertices = [b,a];

% Triangulate all pts in the isocontour and check which trias are in/out

DT = delaunayTriangulation(Vertices);

fc = DT.incenter;

in = interp2(I, fc(:,1), fc(:,2))==1;

% Show the result

figure,imagesc(I), hold on,

patch('vertices',DT.Points,'faces',DT.ConnectivityList(in,:),'FaceColor','g')

patch('vertices',DT.Points,'faces',DT.ConnectivityList(~in,:),'FaceColor','c')

plot(fc(in,1),fc(in,2),'b.', fc(~in,1),fc(~in,2),'y.')

And if you were going for minimal traingulation, you could try something like this:

% Get a reduced set of boundary vertices

bb = bwboundaries(I);

for k = 1:length(bb)

dP = diff(bb{k},[],1);

pdiff = bsxfun(@rdivide, dP, sum(abs(dP),2));

idx = find(any(pdiff - circshift(pdiff,1),2));

bb{k} = bb{k}(idx, :);

end

Vertices = fliplr(cat(1,bb{:}));

% Triangulate all pts in the isocontour and check which trias are in/out

DT = delaunayTriangulation(Vertices);

fc = DT.incenter;

in = interp2(I, fc(:,1), fc(:,2))>0;

figure,imagesc(I), hold on,

patch('vertices',DT.Points,'faces',DT.ConnectivityList(in,:),'FaceColor','g')

patch('vertices',DT.Points,'faces',DT.ConnectivityList(~in,:),'FaceColor','c')

plot(fc(in,1),fc(in,2),'b.', fc(~in,1),fc(~in,2),'y.')

Anand
on 24 Aug 2013

Try using bwperim and delaunay. Something like this:

BW = bwperim(im);

[x,y] = find(BW);

tri = delaunay(x,y);

Hope this helps!

Sathyanarayan Rao
on 10 Aug 2017

Check this code that uses Gmsh

https://nl.mathworks.com/matlabcentral/fileexchange/61507-binary-image-to-finite-element-mesh---gmsh-geo-file--?s_tid=prof_contriblnk

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!