## 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); % for DIPimage 2
% cc = traceobjects(img,31,2,'chain code') % for DIPimage 3, use cc{1} instead of cc.chain below
```

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.