Asked by Geert
on 23 Aug 2013

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')

Answer by Sven
on 28 Aug 2013

Edited by Sven
on 28 Aug 2013

Accepted Answer

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.')

Answer by 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!

Sven
on 25 Aug 2013

I'm afraid this would only work for single shapes that are their own convex hull.

Geert, is your current implementation short enough to post here? This is actually quite a tricky problem, depending on how you want your output, and given your implementation you might get some suggestions of cleaner code to do the same job. There is an "isocontour" file exchange entry that will get a polygon around each of your shapes, however that will not be a set of triangles covering your surface as illustrated in your picture, it will just be the polygon(s) defining the outline.

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.

Answer by 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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.