MATLAB Answers

How to convert binary image to 2D triangulation?

16 views (last 30 days)
Geert
Geert on 23 Aug 2013
Answered: Sathyanarayan Rao on 10 Aug 2017
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')

Accepted Answer

Sven
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.')
  1 Comment
Geert
Geert on 29 Aug 2013
Thx Sven, this provides me with some nice alternative ways for getting the job done. Many thanks! By the way, I improved the robustness of my code as well, I made the proper adjustements in the question...

Sign in to comment.

More Answers (2)

Anand
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!
  2 Comments
Geert
Geert on 28 Aug 2013
Sven,
I've cleaned up my code and put it in a single m-file. I added it to my question. Suggestions for speed-up, improved robustness, etc. are very welcome!

Sign in to comment.


Sathyanarayan Rao
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

Community Treasure Hunt

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

Start Hunting!