Shouldn't bwconvhull() be idempotent?

Shouldn't the following two calculations give the same results?

load BWimage
A=bwconvhull(BW);
B=bwconvhull(bwconvhull(BW));

They don't.

isequal(A,B)

ans = logical
   0

Answer 1

[Steve Eddins]

Hi Matt,

At a general level - and speaking from too much personal experience - expectations about geometrical operations such as the convex hull are often not met when we're dealing with a discrete grid instead of a continuous domain. Since bwconvhull is inherently working on inputs and outputs on a discrete grid, it is really only an approximation of a true convex hull, and I wouldn't expect it to be exactly idempotent.

With approximations, different approximation techniques usually have their own pros and cons, and one can consider whether one approximation technique is better than another, and under what conditions.

Your function treats pixels as points. It computes the convex hull of all the pixel centers and then uses inpolygon to "rasterize" the resulting shape. A disadvantage of that approach is that it will make skinny, single-pixel lines disappear competely because the convex hull of the pixel centers is a degenerate polygon with no area. On the following input, your function produces no foreground pixels at all in the result.

W = zeros(30,30);
W(15,10:20) = 1;

I don't think we could ship this method without somehow addressing that issue.

The function bwconvhull treats pixels as squares having unit area. It computes the convex hull of the set of all the pixel corners and then uses roipoly to rasterize the result. This method behaves better for skinny shapes, but it behaves less well in terms of idempotence, as you have clearly demonstrated.

Oversimplifying a bit, it comes down to tie-breaking rules related to polygons that go exactly through a pixel center. Is such a pixel inside or outside the polygon? Each choice results in its own set of nonideal behaviors.

It is worth asking, though, whether the method in bwconvhull can be improved with respect to idempotence while avoiding the degeneracy I mentioned above. We can take a look at this.

[Matt J]

Thanks, Steve. That does clear things up considerably. However, I wonder if the ideal compromise would be to construct a polygon boundary from the pixel corners (as now), but only count a pixel as inside the hull if all 4 corners lie strictly within the polygon. My code implementation of this below seems to have the same stability behavior as my previous version, but also preserves skinny pixel lines. It is quite a bit slower than bwconvhull() but maybe as an optimized built-in it wouldn't be so bad.

function Ch=bwconvhullMatt(BW)


     d=0.5;
     dI=[-d,-d,+d,+d]; dII=dI*1.0001;
     dJ=[-d,+d,-d,+d]; dJJ=dJ*1.0001;
    persistent I J II JJ
    if isempty(I)
       sz=size(BW);
       [I0,J0]=ndgrid(1:sz(1),1:sz(2));
       I=I0(:)+dI; J=J0(:)+dJ; %list of all pixel corners in image grid
       II=I0(:)+dII; JJ=J0(:)+dJJ; %list of all "dilated" pixel corners
    end

    idx=bwmorph(BW,'remove');
    

    r=II(idx,:); 
    c=JJ(idx,:);
    
    k=convhull(r,c,'Simplify',true);
    r=r(k); c=c(k);
    
    [in,on]=inpolygon(I,J,r, c);
    
    
    
    Ch=reshape( all(in&~on,2) ,  size(BW)  );

end

[Steve Eddins]

Although requiring all four pixel corners to be inside the convex hull does seem to behave better with respect to the idempotence expectation, there are some disadvantages. I give examples below to show that your latest method tends to bias the result about 0.5 pixel in the concave direction along diagonal portions of the convex hull, and that can result in anomalies when processing multiple adjacent objects.

Because of these issues, I don't think it is a clear winner, and I don't think we could justify making an incompatible change to bwconvhull by switching to this method.

A = zeros(30,30);
A(20,5:25) = 1;
A(10:20,5) = 1;
showbig = @(varargin) imshow(varargin{:},'InitialMagnification','fit');
showbig(A)

B = bwconvhullMatt2(A);
showbig(B)
x = [4.5 5.5 25.5 25.5 4.5 4.5];
y = [9.5 9.5 19.5 20.5 20.5 9.5];
hold on
plot(x,y,'g','LineWidth',3)
hold off

You can see a significant average discretization error towards the interior. Compare with bwconvhull:

C = bwconvhull(A);
showbig(C)
hold on
plot(x,y,'g','LineWidth',3)
hold off

The average discretization error is close to zero, and the magnitude of those errors (voids and protrusions) is smaller.

The bias in your latest version can cause the results to not meet other expectations of convex hull when processing multiple, adjacent objects in an image. Here is an example.

W = zeros(30,30);
W(9,5:25) = 1;
W(10:19,25) = 1;
showbig(W)

The objects in image A and image W make a closed box together:

X = 2*A + 3*W;
showbig(X,parula(3))

Your latest method leaves gaps in the union of the convex hulls of the two objects:

Q1 = bwconvhullMatt2(A);
Q2 = bwconvhullMatt2(W);
showbig(Q1 | Q2)

The function roipoly, which is used by bwconvhull, is designed specifically to have better behavior for such scenarios. Compare:

P1 = bwconvhull(A);
P2 = bwconvhull(W);
showbig(P1 | P2)

function Ch=bwconvhullMatt2(BW)
d=0.5;
dI=[-d,-d,+d,+d]; dII=dI*1.0001;
dJ=[-d,+d,-d,+d]; dJJ=dJ*1.0001;

sz=size(BW);
[I0,J0]=ndgrid(1:sz(1),1:sz(2));
I=I0(:)+dI; J=J0(:)+dJ; %list of all pixel corners in image grid
II=I0(:)+dII; JJ=J0(:)+dJJ; %list of all "dilated" pixel corners

idx=bwmorph(BW,'remove');

r=II(idx,:);
c=JJ(idx,:);

k=convhull(r,c,'Simplify',true);
r=r(k); c=c(k);

[in,on]=inpolygon(I,J,r, c);

Ch=reshape( all(in&~on,2) ,  size(BW)  );
end

[Matt J]

Because of these issues, I don't think it is a clear winner, and I don't think we could justify making an incompatible change to bwconvhull by switching to this method.

I don't seriously expect that you would, given the incompatibilities with past bwconvhull behavior.

Note, however, that the flaws you cite in the bwconvhullMatt2 are all based on the idea that the convex region being approximated is the one formed from the corners of the pixels. If instead we compared the result to the convex hull of the pixel centers, then the behavior of bwconvhullMatt2 in your examples is ideal.

I know our criterion for turning a pixel white was based on a larger polygon, but that doesn't necessarily mean that same polygon is the one that the result should be compared against.

[Steve Eddins]

Thanks, Matt. I will take a further look at your suggestions to see if there is something useful we can do.

[Matt J]

Thank you, too. It was a good discussion.

Answer 2

[Jonas]

ideally it should. already the shape you provide in the original is convex, but the result of the first call of bwconvhull and the original differ at 130 positions. the result of the first bwconvhull result and second show differences at 114 positions.

i think the general issue is that shapes/lines can not be represented good in pixel format except for horizontal, vertical and diagonal lines. in the other cases the algorithm just rounds some pixel to be inside the hull or not

using bwconvhull on easier shapes like a block or something like

[0 0 0 0; 
 0 1 1 0;
 0 1 1 1;
 0 0 1 0];

is idempotent

[Matt J]

i think the general issue is that shapes/lines can not be represented good in pixel format except for horizontal, vertical and diagonal lines.

Sure, but the test below shows that whatever method bwconvhull is using, it is strangely and unnecessarily unstable. WIth my own implementation, I am able to iterate the convex hull operation many times without distorting the original shape nearly so much as bwconvhull does.

load tst

B=A;
tic;
for i=1:100; B=bwconvhull(B); end
toc;

Elapsed time is 0.475528 seconds.

C=A;

tic;
for i=1:100; C=bwconvhullMatt(C); end
toc

Elapsed time is 2.404097 seconds.

tiledlayout(1,3,'Padding','none');
nexttile, imshow(A); title 'Original'
nexttile, imshow(B); title 'BWCONVHULL'
nexttile, imshow(C); title 'Matt''s BWCONVHULL'

function Ch=bwconvhullMatt(BW)

persistent I J
if isempty(I)
   sz=size(BW);
   [I,J]=ndgrid(1:sz(1),1:sz(2));
end

 V=cell2mat(bwboundaries(BW));
 warning off
 pgon=convhull(polyshape(V));
 warning on

 Ch=reshape( inpolygon(I(:),J(:),pgon.Vertices(:,1), pgon.Vertices(:,2)), size(I)  );
 
end

[Jonas]

another example where built in matlab functions work worse than user developed functions.

Answer 3

[Matt J]

Below is a 3rd version of the function I proposed in my various conversations with Steve and Jonas. Assuming one is happy with treating the convex hull as the hull of the pixel centers, this version will compute it faster than the previous versions, and without the problems handling skinny 1 pixel lines. It does, however, have the issue with adjacent regions mentioned by Steve, though an imclose() operation could be used to fix that.

load tst

B=A;
tic;
for i=1:100; B=bwconvhull(B); end
toc;

Elapsed time is 0.690813 seconds.

C=A;
tic;
for i=1:100; C=bwconvhullCenters(C); end
toc

Elapsed time is 0.809547 seconds.

tiledlayout(1,3,'Padding','none');
nexttile, imshow(A); title 'Original'
nexttile, imshow(B); title 'bwconvhull'
nexttile, imshow(C); title 'bwconvhullCenters'

function Ch=bwconvhullCenters(BW)



    d=0.0001;
    dr=[-d,-d,+d,+d];
    dc=[-d,+d,-d,+d];

 
    V=bwboundaries(BW);

    [r,c]=deal(V{1}(:,1), V{1}(:,2));
    r=r+dr;
    c=c+dc;
    

    k=convhull(r,c,'Simplify',true);
    r=r(k); c=c(k);


    Ch=roipoly(BW,c,r);
    
end