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.
I = phantom('Modified Shepp-Logan', nx)>0;
temp = zeros(size(I)+2);
temp(2:end-1,2:end-1) = I;
I = temp;
contourThreshold = 0.5;
[Lines,Vertices,Objects] = isocontour(I,contourThreshold);
Vertices = fliplr(Vertices);
DT = delaunayTriangulation(Vertices);
fc = DT.incenter;
in = interp2(I, fc(:,1), fc(:,2))>=contourThreshold;
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:
[a,b] = find(bwperim(I));
Vertices = [b,a];
DT = delaunayTriangulation(Vertices);
fc = DT.incenter;
in = interp2(I, fc(:,1), fc(:,2))==1;
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:
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{:}));
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.')
Best Answer