Cris’ Image Analysis Blog

Last month I wrote a post showing how to calculate the perimeter of an object using its chain code. In this post I want to review several more measures that can be easily obtained from the chain codes: the minimum bounding box; the object’s orientation, maximum length and minimum width; and the object’s area. The bounding box and area are actually easier computed from the binary image, but if one needs to extract the chain code any way (for example to compute the perimeter) then it’s quite efficient to use the chain code to compute these measures, rather than using the full image. To obtain the chain codes, one can use the algorithm described in the previous post, or the DIPimage function dip_imagechaincode [Note: this function no longer exists in DIPimage 3, use the traceobjects function instead.].

The bounding box

The bounding box of an object is easy to obtain: one need only know the minimum and maximum x-coordinate and y-coordinate (note this code requires DIPimage):

img = label(~threshold(readim('cermet')));
binimg = img==33;
coords = findcoord(binimg);
boundingbox = [min(coords),max(coords)]

boundingbox =
   151   114   179   149

Just for fun, let’s draw the bounding box over the image:

img
dipmapping('labels')
hold on
boundingbox(3:4) = boundingbox(3:4)-boundingbox(1:2);
rectangle('position',boundingbox,'edgecolor','red');

Now, if we have the chain code, and the start coordinates for the chain code, we can derive the coordinates of all the boundary pixels. The minimum and maximum of these coordinates are also the minimum and maximum of the coordinates of all the pixels in the image:

cc = dip_imagechaincode(img,2,33);  % DIPimage 2 only
directions = [ 1, 0
               1,-1
               0,-1
              -1,-1
              -1, 0
              -1, 1
               0, 1
               1, 1]; % directions(chaincode+1,:) = [x,y] increment
coords = cumsum([cc.start;directions(cc.chain+1,:)]);
boundingbox = [min(coords),max(coords)]

boundingbox =
   151   114   179   149

Feret diameters

More interesting are the so-called Feret diameters. The name refers back to a publication by L.R. Feret, “La Grosseur des Grains” (Assoc. Intern. Essais Math. 2D, Zurich, 1931). I have never been able to find this publication, though, and only know of it through references in other papers. In any case, it appears Feret defined some size measure for objects. One of them, the width of the smallest projection, yields the smallest possible rectangular box that you can put around an object. The projection perpendicular to it is not necessarily the longest possible projection. Take for example a square object. The smallest and longest projections are 45 degrees apart. To compute the longest and shortest projections, one needs to rotate the object over small angles and measure the projections. The most efficient way to do this happens to be using the chain code. It is much cheaper to compute the coordinates after rotation of only the boundary pixels, compared to rotating the whole object.

There are two ways to compute the coordinates after rotation of the boundary pixels.

1. We take the coords array we computed above, and multiply it with the rotation matrix. This yields the coordinates of the rotated boundary. The maximum and minimum coordinates can now be found, the difference between the two yields the width of the projection.

2. We compute the coords array anew for each rotation, just like above, but using a rotated version of the directions array. This involves fewer multiplications, so potentially is more efficient. This code computes the projection sizes using the second method:

stepsize = 2*(pi/180);
R = [cos(stepsize),-sin(stepsize);sin(stepsize),cos(stepsize)];
maxD = -Inf;
minD = Inf;
for angle=0:2:88
   coords = cumsum(directions(cc.chain+1,:));
   sz = max(coords)-min(coords)+1;
   if maxD<sz(1)
      maxD = sz(1);
      maxA = angle;
   end
   if maxD<sz(2)
      maxD = sz(2);
      maxA = angle+90;
   end
   if minD>sz(1)
      minD = sz(1);
      minP = sz(2);
      minA = angle;
   end
   if minD>sz(2)
      minD = sz(2);
      minP = sz(1);
      minA = angle+90;
   end
   % rotate coordinate system by stepsize
   directions = (R*directions')';
end
fprintf('Object length: %.2f, at %d degreesn',maxD,maxA);
fprintf('Minimum bounding box: %.2f by %.2f, at %d degreesn',...
        minD,minP,minA);

Object length: 38.80, at 128 degrees
Minimum bounding box: 23.05 by 38.29, at 28 degrees

Area

Object area is the most common measure I know of. It is very easy to compute by just visiting all image pixels and counting the number of pixels that belong to the object. On a modern architecture, this is much cheaper than obtaining the chain code. But if you already have the chain code, computing the area from it is very efficient. It is always good not to have to transverse the image a second time. For example, if one is looking for the roundness (variably called “shape factor” and “circularity”, a measure derived from the perimeter and the area), it is possible to trace the object’s boundary once and obtain both the perimeter and area in one go.

To understand the algorithm, imagine a baseline underneath the object. For each boundary pixel on the top of the object, you add the number of pixels in between it and the baseline. For each boundary pixel on the bottom of the object, you subtract the number of pixels in between it and the baseline. The net will be the number of pixels between the top and bottom boundary. The beautiful thing is, that it doesn’t matter where this baseline is, it can be at any height, it can even cut the object or be above the object. The math stays the same. As you move from one boundary pixel to the next, you evaluate whether you are on the top or bottom of the object (if you’re going to the right, you’re on the top), and add or subtract the current value for the vertical distance to the baseline. At the same time, you adjust this distance as you move up or down in the image.

This algorithm is due to P. Zamperoni (Signal Processing 3(3):267-271, 1981), though due to a rather unspecific paper I had to derive the exact increments anew. Download the function, and try to understand how it works: ccarea.m. The value B is initialized randomly. This is the vertical distance between the first pixel in the chain and the baseline. I initialized it here to 1000 just to make a point, but you could choose 0 just as well. Within the loop, the algorithm first corrects for special cases (these are given by the specific sequence of two consecutive chain codes), then performs the “normal” part of the algorithm (which depends only on the current chain code): a step to the right adds to the area measure, A, and a step to the left subtracts from it; a step up increases the distance to the baseline, B, and a step down decreases it; diagonal steps are a combination of a vertical and horizontal step, so affect both variables. To understand the correction part of the loop, think of the pixels in the leftmost column of the object: these are never added to the calculated area. Similarly, the pixels below the rightmost column are never subtracted. The correction simply takes care of these “lost pixels” as the chain code changes direction.

As you can see, the algorithm produces the same result as sum on the object:

ccarea(cc.chain)
sum(binimg)

ans =
   641

ans =
   641

If the object has a hole in it, the chain code will not be able to represent the object correctly. The area computed from the chain code will be that of the object with all holes filled. This will, obviously, be different from the result of sum.