## The convex hull of a 2D object

Last year I wrote about computing the the boundary length and various other measures, given an object’s chain code. The chain code is a simple way of encoding the polygon that represents a 2D object. It is very simple to compute the object’s convex hull given this polygon. Why would I want to do that? Well, the convex hull gives several interesting object properties, such as the convexity (object area divided by convex hull area). Certain other properties, such as the Feret diameters, are identical for an object and its convex hull, and the convex hull thus gives an efficient algorithm to compute these properties.

Avraham Melkman (1987) has published a very efficient algorithm to compute the convex hull given a polygon (find the PDF through Google Scholar). It has been said that this might very well be the most efficient algorithm possible. It visits each polygon vertex only once, and keeps only those vertices that form the convex hull. To accomplish this, it uses a deque (also known as double-ended queue) to represent the convex hull. The deque has 4 basic operations: *push* and *pop* add and remove elements from the end of the queue, and *insert* and *remove* add and remove elements from the front of the queue. By using such a data structure, Melkman was able to add and remove vertices from his convex polygon at either side. The points in the deque always constitute a convex polygon, and is the convex hull of the input points seen up to that point.

The Melkman algorithm also depends on a function, which I’ve called *position* here, which returns 1 if `pt3` is to the right of the directed line formed by `pt1` and `pt2`, -1 if it is to the left, or 0 if all three points are collinear (form a straight line). This function simply returns the sign of the signed area of the parallelogram given by the vectors `pt2-pt1` and `pt3-pt1`. It we embed these vectors in 3D (adding the value 0 as the z component), and compute the cross product, the norm of this cross product gives the area of the parallelogram, whereas its direction is given by the ordering of the three points on the plane (as per the right-hand rule). The cross product will, however, always have x and y components equal to 0, so we only need to compute the z component. Let’s define an anonymous function `position(pt1,pt2,pt3)` for this:

position = @(pt1,pt2,pt3) ... sign( (pt2(1)-pt1(1))*(pt3(2)-pt1(2)) - ... (pt2(2)-pt1(2))*(pt3(1)-pt1(1)) );

Before we can apply the algorithm, we’ll need some data. I like the image “cermet”. Let’s pick object #31, and get its chain code:

```
img = label(readim('cermet')<120);
cc = dip_imagechaincode(img,2,31);
```

The code below, however, takes a polygon represented by a set of points, not a chain code. We can easily convert the chain code `cc` into a list of coordinates `coords` (this is the same code as in the post on chain codes):

directions = [ 1, 0 1,-1 0,-1 -1,-1 -1, 0 -1, 1 0, 1 1, 1]; coords = cumsum(directions(cc.chain+1,:)); N = size(coords,1); % Draw input polygon: figure('position',[1000,500,450,450]); plot([coords(:,1);coords(1,1)],[coords(:,2);coords(1,2)],'b-') hold on, axis equal, axis off

Melkman’s algorithm starts by taking the first three points in the list and forming them into a start of the convex hull. However, the algorithm assumes these points are not collinear. So let’s examine the first 3 points in the set, and delete the second one if all 3 are collinear:

while N>=3 && position(coords(1,:),coords(2,:),coords(3,:)) == 0 coords(2,:) = []; N = size(coords,1); end

At this point, we’d of course check to make sure that we have more than 3 points left. If not, we’re done! There is also a special case that could happen if our object is a single line. In this case, all points in the object are collinear, and the code above leaves only 2 points. We’d need to detect this case, and make sure that we keep the two extreme points, rather than the first and last points in the list. Otherwise length measurements on our convex hull will be wrong.

We start by creating a deque data structure `D` (in MATLAB we can just use a simple array), and put the first three points in it. We first check the order of these points, and put them in the deque so that the polygon is right-turning. Furthermore, we always want the first and last point in our deque to be the same point:

if position(coords(1,:),coords(2,:),coords(3,:)) > 0 D = coords([3,1,2,3],:); else D = coords([3,2,1,3],:); end % Add convex hull to figure: h = plot(D(:,1),D(:,2),'ro-');

Next, the algorithm takes one new point `v` at the time, and checks where it lies with respect to the first two points and the last two points in the deque. This tells us if the point is already inside the convex hull or not. If it is inside, we ignore the point and take the next one. If it is not inside, we’ll be wanting to add this point at the beginning and at the end of the deque (push and insert). But before we add this new point, we need to remove (pop and remove) previous points that will potentially end up inside the convex hull. We start with the end of the deque, and see if the last point is inside the line formed by the second to last point and the new one. If so, we pop that point. This is repeated until we cannot pop the last point any more, and then push the new point. We repeat the same procedure with the front of the deque, removing points and then inserting the new one. Now the algorithm loops back to process the next point. This is the code:

ii = 3; while 1 ii = ii+1; if ii>N,break,end v = coords(ii,:); while position(v,D(1,:),D(2,:)) >= 0 && ... position(D(end-1,:),D(end,:),v) >= 0 ii = ii+1; if ii>N,break,end v = coords(ii,:); end if ii>N,break,end while position(D(end-1,:),D(end,:),v) <= 0 D(end,:) = []; % pop end D = [D;v]; % push v while position(v,D(1,:),D(2,:)) <= 0 D(1,:) = []; % remove end D = [v;D]; % insert v % Update convex hull in figure: set(h,'xdata',D(:,1),'ydata',D(:,2)); pause(0.2) end D(end,:) = []; % start and end points are the same, remove one.

When running the code, you will see the convex hull being updated as new points are examined by the algorithm (note that the slowness comes from the `pause` statement in the code!). At the end of the algorithm you’ll see this on your screen:

In a future post we’ll use this convex hull to compute Feret diameters.

Dear Chris,

thank you for the many useful tricks that you publish here!

My idea is to use convexity and chaincode to find locally concave and convex parts of my boundary. My objects can look similar to your example boundary. More specifically I want to look at circular-like objects that touch and morph together giving two circles with different degrees of overlap/pinching. I am hoping that I can identify regions where circular parts connect, because this is where the local curvature becomes concave. Then I want to split my object’s boundary into the circular sub-parts composing the individual circle-like parts.

Similar to your chaincode here, my boundaries are not smooth enough/too noisy to look at “diff(cc)” directly (additionally I am not sure how to find the curvature when the direction jumps from 7->0).

Maybe one can identify the concave regions by looking at parts of the boundary that are far away from the boundary of the convex hull? The actual concave part of the boundary can be very short as well, such that a plot of “cc” shows a noisy, monotonously decrasing function, with discontinuities at regions where the boundary goes from one circle to the next circle.

Do you have any comment/idea on how to measure local curvature of a boundary and or how to find discontinuities in the local average curvature?

Many thanks for any advice on this!

Hi Enno,

This is an old problem, much code has been written to try to solve it. One thing you could do is take the coordinates of the boundary pixels, in order (as you’d obtain from the chain code), and smooth the two vectors formed by the x-coordinates and the y-coordinates. That is, given 20 pixels along the boundary, you’d have two vectors with 20 values each. You can smooth these by applying a low-pass filter (preferably a Gaussian), and using periodic boundary condition (i.e. assuming

`x(end+1)=x(1)`

). From these sequence of points you can compute the second derivative (extra points if you do the smoothing and derivative with the same filter!, see an earlier post on Gaussian filtering). The second derivative will be negative everywhere except at two points, hopefully.However, it should be much easier to separate your circles using the distance transform + watershed trick. From your binary image, compute the distance transform; this gives the distance to the background for every object pixel. Two circular objects, even when somewhat overlapping each other, will each have a local maximum in their middle. The inverse of this distance transform can be used as input to the watershed, which will create one region for each local minimum (hence using the inverse, so that each object has a local minimum in its middle). The watershed will create a boundary that is exactly like the one you are looking for.

Dear chris

i have some binary images that want to cllasify them base on there shape .if they have circular or elliptical shape they belong to class one,if they be ellipse with dent in their boundary they belong to class one. i dont know how can use this feature. can you help me how can i do this?

Thank your for your helping in advande…

Zeinab,

You could probably compare the object’s main axes lengths (computing perimeter of an ellipse with that size) with its measured perimeter. Axes lengths measured using Feret diameters, perimeter measured using chain code algorithm. Both methods are explained elsewhere in this blog.