# How to convert binary image to 2D triangulation?

조회 수: 15 (최근 30일)
Geert 2013년 8월 23일
답변: Sathyanarayan Rao 2017년 8월 10일
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 2013년 8월 28일
편집: Sven 2013년 8월 28일
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이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Geert 2013년 8월 29일
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...

댓글을 달려면 로그인하십시오.

### 추가 답변 (2개)

Anand 2013년 8월 24일
Try using bwperim and delaunay. Something like this:
BW = bwperim(im);
[x,y] = find(BW);
tri = delaunay(x,y);
Hope this helps!
##### 댓글 수: 2없음 표시없음 숨기기
Sven 2013년 8월 25일
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 2013년 8월 28일
편집: Geert 2013년 8월 28일
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!

댓글을 달려면 로그인하십시오.

Sathyanarayan Rao 2017년 8월 10일
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

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Object Analysis에 대해 자세히 알아보기

### Community Treasure Hunt

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

Start Hunting!

Translated by